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ABSTRACT 

Projection algorithm has "been found to have a number 
of desirable properties which make the algorithm attractive 
for practical applications. In this work the algorithm has 
been reviewed in the context of reconstruction problem. 
Effects of binary orthogonalization on the speed 'of con- 
vergence of this algorithm have been numerically studied 

Vp,' 

using a set of 32x32 digital functions. It is shown that 
orthogonalization does expedite convergence of the pro- 
jection method and results in a securing of about 30-50$ 
in computational time. 



CHAPTER 1 


INTRODUCTION 

The 1979 Nohel Prize for Medicine and Physiology was 
shared by Allan M Cormack and G N Hounsfield for the 
development of Computer A ss i s ted Tomography (CAT), a new 
method for X ray diagnosis. This non invasive diagnostic 
procedure makes use of computers in imaging internal organs 
within human body by mathematically combining X-ray image 
data from numerous angles to obtain the image of any cross- 
section within the body. The technique involved is ’the 
reconstruction from projections’. This technique has been 
applied, in various forms, to other fields of scientific 
endeavour such as radio astronomy, microscopy, biomedical 
engineering, cellular biology, nuclear science, digital 
image transmission etc. The list is limited only by imagi- 
nation. 

In most of the cases images are used to represent 
the distribution of some property of an object or physical 
system e.g. molecular structure, X ray. Most images are 
formed directly by optical instruments using visible 
light reflected or transmitted by an object. In many 
applications in which an image is required, only indirect 
measurements by probing an object can be made with visible 
radiation or by interpreting radiations emitted by it. 
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Often the measurement data is not in a suitable form for 
immediate interpretation. The general aim of all image 
re cons tract ion procedures is to process the data to form an 
image so as to facilitate the interpretation of the mea- 
surement . 

A number of techniques for reconstruction from 
projections are available and a brief survey is given in 
the next section. 

1.1 A Word on Present Reconstruction Techniques 

J. Radon did the pioneering work in the field of re- 
construction and proved that it is possible to reconstruct 
structures from their projections. The algorithms which 
have come up in the last few years can be classified into 
the following : 

1. Direct matrix inversion techniques [31(421 

2 . Summation, linear superposition, back projection [4] [4 3 J 

3. Algebraic reconstruction techniques (ART) (51 [44 J 

4 . Simultaneous iterative reconstruction technique (SIRE) [6J 

5 . Orthogonal tangent correction E 71 

6. Iterative least squares techniques [3] [ 8 ] 

7. Convolution techniques f 3][9J [451 
S. Geometric mean iterative techniqaes 
9. Pourier reconstruction dll 1 46J 

10. Summation of projections after Hilbert transform of 

the derivative of the projections. 1221 [481 [491 [501 
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The techniques mentioned above can be grouped broadly into 
the following four categories : 

1. Summation (1) 

2. Use of Fourier transforms (7,9) 

3. Analytic solution of integral equation (7,9,10) 

4. Series expansion approaches (3, 4, 5, 6) 

(The numbers in small brackets refer to the serial number of 
the technique given above ). 

More generally, these techniques may be grouped as : 
signal space reconstruction and Fourier space reconstruction. 

Bracewell and Fiddle [141 proposed algorithm for recon- 
struction in signal space while De Rosier and A. Klug [ 11] 
did the initial work in the area of Fourier transform 
methods. The Direct Matrix Inversion techniques were dealt 
with by Budinger and Gullberg [3] [13] who also showed that 
this method requires large storage and long computational 
times. Andrews discussed summation, linear superposition, 
back projection 1 41 method. Then Vainshtein [16 J 
Ramachandran and Lakshminarayan [20] and Bracewell and 
Riddle f 9l came forward with the convolution method. 

Gordon, Bender and Herman [51 and Gordon [21] proposed 
Algebraic Reconstruction Technique which they claimed to be 
the best (under certain circumstances) among the contemporaiy 
techniques. It was implemented in signal space. Gilbert 16 ] 
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came forward with iterative relaxation methods. Krishnan, 
Prabhu and Krishnamurthy [ 15 ] gave an iterative scheme 
"based on unbiased probabilistic estimates. Geometric mean 
iterative technique was given by P. Schmidlin [10] . D.E. Kuhl 
proposed orthogonal tangent correction method. In 1971 
Yainstein came up with a summation method. Budinger and 
Goullberg [131 implemented simultaneous reconstruction 
technique (SIKF). Wernecke and D'Addario H7l have applied 
the principles of maximum entropy to the two dimensional 
digital image reconstruction problem. 

1.1.1 Choice of Algorithm for our Study 

As we have seen above a large number of algorithms 
exist for reconstructing three, dimensional or two dimen- 
sional images from their projections. We needed for our 
study an algorithm which is comparatively well behaved, iie. 
its computationally less cumbersome, essentially performs 
pseudoinversion and is guaranteed to converge. Not much 
has been said about the relative merits and demerits of 
various algorithms and we base our decision on whatever 
comparison is available. 

Bellman, Bender, Gordon and Rowe [371 have shown that 
ARP does comparatively better than the Pourier techniques 
particularly when a limited number of projections in a 
limited range of angles, is available. Herman [38] compared 
an improved summation method with ART for 4 test patterns, 
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4 distinct sets of projections and 18 different criteria 
of success. ART -was consistently better. Herman and 
Rowland [391 compared ART with SIRT and convolution method. 

It was found that for und er determ in ed , noisy system for 
which projections are available from a narrow range of angles, 
ART is superior. Schmidlin f 40l found ART to be superior 
to the summation method. Gordon compared convolution and 
a form of ART and the latter was found to be far superior [411. 

The Projection Algorithm (PRA) which was first proposed 
by Kacmarz [30] as a method of solving linear equations and 
later reviewed by Tanabe (311 who proved its guaranteed 
convergence is comput ationally , in principle, the forerunner 
of ART. Recently it was applied to the problem of restora- 
tion by Huang and Ramakrishna [331 and was found to be very 
well behaved. 

Keeping the above considerations in mind it was felt 
that PRA could profitably be utilised for the reconstruction 
problem. Accordingly, this thesis is concerned with compu- 
tational study of the PRA for the problem in hand. 

1.2 Present Work 

The reconstruction algorithm which we have mentioned 
above reconstructs images from their projections with 
varying speed and accuracy. Some of these algorithms are 
slow and are suitable only for batch mode processing while 
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others are faster and allow real time processing. Even at 
present most of the algorithms are either Quite slow or 
else they require prohibitively large storage. Many pra- 
ctical applications require faster algorithms for the 
result to be useful. 

Ramakrishna 133] advocated the use of some form of 
orthogonalization with the projection method to expedite 
its convergence. In the context of the restoration problem 
it was shown comput ationally that binazy orthogonalization 
followed by PRA resulted in considerable savings in compu- 
tational time. Similar techniques have been studied, in 
this thesis, in the context of reconstruction problem. PRA 
with orthogonalization, or simply OPRA, has been compared 
with the PRA to show formers superiority in many respects. 
The bases of comparison are the number of mismatches, number 
of projections, number of iterations and some other per- 
formance criteria. It has been shown that the gain in 
computational speed are in the range of 20$ - 50$. 

1 .3 Overview of the Thesis 

Chapters 2 and 3 deal with the basic ideas of recon- 
struction problem in two dimensions. They discuss various 
steps involved in solving the reconstruction problem digi- 
tally. These include digitization, obtaining projections 
and the formulation of the problem as a set of a large 
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number of simultaneous linear equations. Chapter 3, in 
particular, presents the problem in the form in which it will 
be taken for solution in the later chapters. 

Chapter 4 describes the projection method and its speed 
up by an orthogonalization scheme. Binary orthogonalization 
has been discussed in some detail. 

Chapter 5 contains the results obtained for the test 
pictures used. A comparison between PRA and OPRA has been 
made to demonstrate the superiority of OPRA over PRA. 

Appendix A' discusses reduction of dimensionality of the 
reconstruction problem, reduction of the storage requirement 
and of takes up the question of compactness of figures in 
digital pictures. 

Appendix B includes a discussion on the number of views 
required for reconstruction. It also deals with the limits 
on discretization of projection angles. 

Appendix C gives the listing of the combination of 
characters used to obtain grey levels in time printer imagery 
The list of the Fortran programs used in computations is 
contained in Appendix D. 
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CHAPTER 2 

ON RECONSTRUCTION FROM PROJECTIONS 
2.1 Introduction 


J. Radon, an austrian mathematician, triggered the 
development in the field of reconstruction by his historical 
assertion made in 1917 that a two or three dimensional structure 
could be reconstructed uniquely from the infinite set of all its 
projections. Since then, this result has been exploited in 
developing algorithms for reconstructing images from its proje- 
ctions. Many of the algorithms are mechanisms of evaluating 
Radon's relation between the value of each picture element in 
polar coordinates A(r,©) and the projections for all angles 
fe(x,0), Fig. (2.5), where x denotes an element along the pro- 
jection corresponding to the line integral through the section 
to be reconstructed. Thus Radon [22] in 1 91 7 and subsequently 
others [12] t [49] showed that 

»/2 


A(r,0) 


1 

m * . ,m wpmi 

9 TT 




JT/2 - 


X 


“sin (9-©' ) — x dx d0 ’ ^ 2>1 ^ 


That this method Is practically non-implementable can be seen 
by the requirement of projections at 'all' angles, i.e. infinite 
number of projections. In practical situations there are only 
a finite number of views and each measurement is subject to 
errors. Over the last twenty years specialized techniques have 
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been developed for solving the problem of estimating the 
distribution of some property in three dimensional space from 
many views or projections at various discrete angles. But not 
many techniques allow an effective compromise between compu- 
tational time and storage requirements. 

2.2 The Versatile Tool 

Computer Assisted Tomography (CAT), which won Dr&. Allan M 
Cormack and GET Hounsfield their Nobel prize, involves recon- 
struction from projections. The technique of reconstruction from 
projections has been applied, in various forms to other fields 
such as computer science, theoretical molecular and cellular 
biology, medical and biological engineering, radiology, 
optics, crystallography , physics, nuclear science, material 
science, electrical engineering to name a few. 

Algorithms which reconstruct three dimensional object from 
two dimensional slices are of great importance in biology for 
they allow us to discover the structure of macromolecules and 
macromolecular assemblies such as ribosomes [ 23 , 24 ] . In 
astronomy we can reconstruct the 2-D images of the galaxies from 
radio signals reaching the earth [ 141 . In the field of medicine 
and physiology this method finds a very important application, 
viz. ^detection of tumors and similar growths in human body. 
Projections are taken at various angles and an image constructed 
from these brings out the tumor as an area of different (greater) 
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density as compared to the surrounding tissue. A three 
dimensional reconstruct can give the size of the tumor and its 
location with amazing accuracy. Patient studies for isotope 
distribution in the head, heart and liver can be accomplished 
by rotating the patient before a scintillation camera in 10° 
increments. The study time is approximately 30 minutes and the 
doses are no greater than the routine studies with 100 to 400 
mrads. Application of these techniques to the heart involves 
gating the camera or the computer to overcome motion M 3 J. 

Another important application is in the area of digital 
image transmission. Application of the concept of projections 
results in a considerable saving in both hardware and proce- 
ssing time. 

2.3 Basic Steps in Reconstruction from Projections 

We will study here the basic theory behind all techniques 
of reconstruction using projections before going on to a parti- 
cular method in later chapters. In the context of this general 
discussion it is sufficient to define a picture as a real valued, 
non negative function of two variables. The value of this 
function at any point gives the grey level at that point [28 3. 

We are considering pictures with grey levels only. 

A projection has been defined as an ordered set of values, 
each being the integral or sum of grey levels along one of a set 
of parallel rays which are bands of finite width drawn across the 
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picture (Fig. 2.1). The set { b^bg, ... b^} of the ray sums 
forms a projection. We will talk about projections and method 
of obtaining them in the next chapter. 

The problem is to reconstruct from a finite subset of its 
projections taken at distinct angles, the picture under test. 

The amount of data available is usually less than that required 
for an exact solution and the use of iterative techniques is 
preferable 126 1 if digital processing is to be employed. 

2.3.1 Discretization 

For digital reconstruction we must discretize the problem. 
In other words we must sample our picture. One way to sample 
it is to break it down into small square cells of equal di- 
mensions with the help of a grid of appropriate size. This 
method is discussed in greater detail later. Another way is 
to sample the picture with an array of equally spaced points 
(Fig.. 2.2). Any method of sampling will yield some discrete 
values which can be expressed in the form of an array. 

The difficulties faced in practical digitization will be 
discussed in the next chapter. 

With the help of Fig. 2.2 we can understand how, after 
digitization, the samples are obtained and how the projections 
are taken. The figure shows an irregular outline which is the 
two dimensional cross-section of our three dimensional density 
distribution or activity. The greater the density at any point, 
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greater is the activity at that point. The figure also shows, 

on the plane of cross-section, a nxn array of sampling points 

which are used to discretize this analog density distribution. 

The density distribution of the object determines the activity 

2 

to be assigned to each of the n points of the sampling array. 

The point falling in the region of greater density gets higher 
activity attached to it and so on. We will see in the next 
section how we use this discretized distribution for obtaining 
projections. The terms density and activity will be used synony- 
mously in the text. 

2.3.2 Obtaining Projections 

To get projections, the activity in the plane of cross- 
section i.e. the activity associated with all the points, is 
projected onto a line (OP in Pig. 2.2) at angle 0^ to some 
reference line (OR in Pig. 2.2). One can see in the same figure 
a set of parallel rays <5^, ^ which are all parallel to 
each other and to the plane of cross-section. Wot much obser- 
vation is required to make out that any of the sampling point 
can fall, in general, in only one of these rays. If it falls 
on the edge then it is associated with the ray to which it 
contributes most. The sum of activities of all the sample 
points in any ray is known as the ray sum of that ray. To be 
more specific we can think of a projection as being formed by 
a group of parallel rays. The projection value associated with 
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ray kO is bg^. which is the sum of the activities of all the 

points included in the ray. The set of these values form a 

projection. Other projections are got by projecting the 

activities a. . onto the line OP at different angles to the 
i j 

same reference line OR. The set bg^, i = 1, ... n forms a 
projection at an angle 0. 


The problem now boils down to reconstruction of the true 

values of the activities a. . (i:1,n, j:1,n) from these pro- 

i j 

jections. We will see in the next few sections how the problem 
can further be reduced to a problem of solving a set of simul- 
taneous linear equations and later, in the next chapter, we will 
formulate the problem in the manner relevent to our discussion. 

2.3.3 Iterative Schemes 


In the iterative methods a. - are continually reconstructed 
in the hope of obtaining a better approximation of the true 
values on each successive iteration. Some algorithm is developed 
to calculate the reconstructs a^ from bg^. We start with 
assumed values of the activities a. . and update these values 
using the projection values bg^. The values of a^ obtained 
above are projected again onto lines at different angles to 
produce a set of reconstructed projection values bg^’s. a ij' s 
are then recomputed according to some function of bg^ and bg^ 
to give the values for the next iteration. The process is 
terminated when satisfactory reconstruction has been obtained. 
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2.4 The Reconstruction Problem as a Set of Simultaneous Equation 

The reconstruction space is represented by a basis set which 
subdivides it into nxn non overlapping subregions of similar 
nature. The unknown density distribution is approximated by the 
values assigned to each element by the reconstruction algorithm. 

Each projection which is finite in extent is also divided 
into non overlapping elements the maximum size of which is 
dictated by the presumed spatial resolution in the projection. 

A projection consists of bands of parallel radiation 
through the density distribution at a known angle 9^.. Each band 
is known as a ray and covers some region of the distribution. 
The projected elements are b.'s (the information about the 

J 

angle has been omitted for ease in writing and bg^ is written 

simply as b,) and corresponding to each b. there is a 5 . 

J J J 

depending on the path of the radiation (Pig. 2.1). b. is a 

J 

function of the sum of the activities of all the points lying 
in the corresponding region <$. . The projected elements 

v 

b-, j = 1, ... n for any 9, form a projection at that angle. 

J K 

If we take projections at p different angles (9^., k = 1, ... p) 

then we have b.., j = 1, .... np projection elements in all, 

J 

each having a corresponding . 

J 

If we denote our reconstruct ion space by S and x(r) as 

our unknown density at point r in S then the integration of 

all such densities in any region 6- will yield approximately 

J 
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the measured value of the projection element b.. Mathemati- 

J 

cally speaking 

j x(r) dr = b j j = 1 , . . . np (2.2) 

j 

The above equation is the starting point of all reconstruction 
methods using projections. 

If -we denote s^ as the ith element of S then -we have 
2 

s*, i = 1, ... n elements. Any section 6. will contain a 

J 

number of such regions (or parts thereof). Then the integral 

in (2.2) over 6. can be broken down into summation of a number 

J 

of integrals over the regions s. which lie in <5.. Thus, 

J 

\ x(r) dr = bj can be written as 

r 

•3 

2 

n 

x(r) dr = b.j j = 1, ... np 

^ s i * i = 1, ... n 2 (2.3) 

This can be arrived at by the definition of a projection 

element, we gave earlier, as the sum of all the densities- 

associated with the corresponding ray. 6. f\ s- isolates 

J 

those cells which contribute to the ray intercepting region 

6.= of S and having the projection element b .. 

J D 

To start with, the density x(t) is unknown and any algorithm 

designed for reconstructing the density distribution will aim 

at approximating the unknown density x(r) by assigning an 
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estimate x^ of its value to each region Sj_. The best 
estimate would result when x^ is the average value of x(r) 
over the subregion Sj_. 


x i = 


[ Si x(r) dr 

? _ dr 

- s i 


i = 1 , . . . n* 


(2,4) 


This ideal outcome cannot be achieved due to limitations on 
the amount and quality of data as well as the reconstruction 
algorithm themselves. We can, nevertheless, get so close to 
this estimate, after making some practical assumptions discu- 
ssed later, that the reconstructed density distribution comes 
very much within the limits of acceptability. 


We do not know how the function x(r) varies within any 

region s^. Without this knowledge we cannot predict the value 

of x(r) in 6^0 s^, the region of the subregion s^ contributing 

to the projection element b - . To go around this difficulty 

J 

we assume that x^ is the average value of x(r) over s^. We 
define a geometric function as 


a . . = I dr i = 1 , ... n‘ 

13 Ji.ns. 

3 i 


(2.5) 


as the area contributing to the element b.. . 
in (2.2) can be approximated as 

x(r) dr 


Then the integral 


I 


s 3 a s -: 


a. . x. 
ID i 


( 2 . 6 ) 
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The equation (2.3) reduces to a set of simultaneous equations : 

2 

n 

I x ± = b.., j = 1, np (2.7) 

i=1 . . 2 

i=1, . . . n 

The above system of equation has the following peculiar 
oharact eristics. 

(1) All the elements of matrix {a. .} are nonnegative. 

X J 

(2) The size of (a. .} is enormous e.g. for n = 512, p = 10 

-A- J 

(a. . )is of the order of 10^. 

X J 

(3) {a. .} is quite sparse since 6.0 s . = 0 f 0 r most pairs (i,j). 

x J J 

(4) The rank of the matrix {a. -> is unknown. 

x J 

(5) The unknown function x(r) is ordinarily assumed to be 
nonnegative so that one desires a solution for which x^ - 0. 

(6) The errors in the data may cause the equations to be 
inconsistent . 

(7) The approximations by which we arrive at (2.7) introduce 
systematic errors which have yet to be analyzed. 

The various questions which might arise at this stage are 

(1) Is it always possible to reconstruct in the above manner? 

(2) Is there any limitation on the number of projections 
required? (3) Is there any limitation on the range of angle 0? 
The discussion below attempts to answer these questions. 
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In some of the following lines we shall prove [ 51 
that if the reconstruction space containing a finite object 
is represented by either nxn squares or nxn points, then n 
projections within an arbitrarily small range of angles suffice 
to reconstruct the space exactly provided only that these are 
also known exactly. 

SEfie fig. 2.3 shows a grid of nxn sampling points used for 
reconstruction of some distribution. We start with defining a 
ray as a band of finite width across the reconstruction array. 
All the rays of a projection are parallel to each other, as 
mentioned earlier. 

The points are represented by an ordered pair (i,j) where 
i denotes the column number and j the row number. 

The first ray of the first projection is chosen such that 
its lower edge passes through a point f^ on the line through 
the second row and also between the points (1,1) and (2,1) (see 
fig. 2.3). This ray now has all the information about the point 
(1,1) and can thus be used to determine the activity associated 
with this point a^,-j • Now choose the first ray of the 
second projection in a similar manner and constrain its lower 
edge to pass through fj and also between the points (2*1) and 
(3^1). This ray will now include the information about the 
activities a^ ^ and a 2 -| at (1,1) and (2,1) respectively. But 
the activity at (1,1) has already been determined so we can 
readily determine the activity a£ -| associated with the point 
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(2,1). By induction we may conclude that n rays (all the 

first rays of n projections) so defined will uniquely determine 

a 1 through a i . The second ray of each projection is 

similarly defined except that it passes through the point F 2 « 

The rays will intercept points of the first row along with the 

points of the second row. But the first row has been completely 

determined so the activities associated with the points in the 

second row can also be determined. In this manner we can 

define n points F . , j = 1, ... n which will give us n rays for 
i J 

each of the n projections, these n projections comprising 
n rays each suffice to determine an nxn picture uniquely. 

Another way of discretizing an analog density distribution 
is by sampling it with a square grid. In this case the density 
distribution is broken up into nxn square cells of same size 
(Fig. 2.4). Each square, in this case, is assumed to have a 
uniform density inside. In this the points F^ , F2, . . . F fl are 
chosen as shown. The cells are denoted in space by the set (i,j) 
where i denotes the column and j denotes the row number. 

The first ray of the first projection passes through the 
point F 1 and top right hand corner of the cell (1,1). This 
ray in doing so intercepts a known fraction of the square (1,1). 
From this information the density associated with this cell 
can be determined. The first ray of the second projection 
passes through F-j and the top right hand corner of the square 
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(2.1) . This ray intercepts a known fraction of the square 

(1.1) and a known fraction of the square (2,1). Then density 
of square (1,1) has already been determined and by the infor- 
mation provided by this ray we can determine the density asso- 
ciated with the square (2,1). Similarly we can define first ray 
of the projections 3,4, ... n which will successively give the 
densities associated with the squares (3,1), (4,1) ... (n,1). 

Now point Pg can taken and from this we can pass the second 
rays of all the projections which also pass through the top 
right hand corner of squares (1,2), (2,2), ... (n,2). This 
determines the activities of the cells (1,2), (2,2), ... (n,2) 
respectively. Similarly we can take points ... P n and get 

n rays of n projections which will together determine the 
activities associated with all of the nxn cells of the array. 

2 

In our study we have assumed the sampling by a grid of n 
cells such that the density associated with each of the cell is 
concentrated at its centre instead of being uniformly distributed 
all over it. This makes the formulation of the reconstruction 
algorithm simpler. We still have a unique solution if n rays 
of n projections are known. 

2.5 Limits on the coarseness of Division of S-the Factor 
Determining the Size of (a. 

tJ 

This problem is taken up in a slightly different context 
in the next chapter where we show how the size of the sampling 
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grid can affect the picture. With coarser grid a circle is 
shown to have been converted into a square or a cross on 
discretization. As the fineness of the grid is increased the 
picture acquires finer details and more details will he 
correspondingly retained in the reconstruction. 


On the basis of the above argument one is convinced that the 
division can be made infinitely fine to get the exact reconstru- 
ction. But since the size of the matrix {a. .} depends directly 
on the number of subregions Sj_ in S, the division should be as 
coarse as possible to restrict the size of s to make computations 
easier. Again -we cannot go on making the division coarser since 

the size of the projection element b - puts a limit on it. How 

d 

bjj restricts the degree of coarseness can be seen from the 
following argument. 


Let c^ be the centroid of the subregion s^ such that c^ e 
(convexity). We define a matrix (u. .} such that 

X d 

1 if c. G 6 g. 

i g a- 

= (2.S) 

0 if not 


The product sum l a, . x. gave the sum of the densities 

X d X 

assigned to all the elements in the region <$.. We can approxi- 

d 

mate the product a. . x. by u. . x. so that 

X d X X d X 


l 

i=1 


U ij X i 



1 , 2 ... 


np 


( 2 . 9 ) 
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Summing up what we have done Jill now since \.e started with 

the equation (2.2) t we have, step by step, approximated this 

integral by the sum of the estimates at a finite set of points 

within 6.. For a given passage <5 the number of points it 

J 

intercepts 

2 

n 

IT. = I u . 3 = 1 , . . . np (2.10) 

O J- J 

i=1 

must be nonzero. This limits the coarseness of the division 

of S. 



P 
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CHAPTER 5 

DIGITAL RECONSTRUCTION - FROM PROJECTIONS 
3.1 Introduction 

In the second chapter we developed the theoretical basis 
for our problem of reconstruction. We discussed as to how 
the projections can be taken and used for reconstructing 
images in two dimensions. It was mentioned that some sort of 
sampling is to be done. In this chapter we look into digi- 
tization more critically outlining the method we are using, 
viz., sampling with a square grid. We also show how pictures 
can be represented in the form of matrices. Once the ques- 
tion of digitization has been settled we develop the problem 
on hand, in the form in which it is going to be solved. 

For convenience we define our analog picture again as 
being a real valued, nonnegative function, f(x,y), of two 
real variables x and y. The value of this function at a 
point is the grey level of the picture at that point (28). 

It is further assumed that this function is well behaved, 
i.e. integrable, Fourier transformable etc. When colour is 
involved the function should be considered as vector valued 
or several functions should be used. We shall confine our 
study to pictures with varying intensities of grey level 
only. 
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We now look into digitization of this analog picture 
so that we can express our problem in a convenient form 

3.2 The Digital Picture 

To achieve representation of continuous space in a 
discrete form many truncated basis sets have been used to 
represent the higher dimensional space in reconstruction 
problem [231 . 

In our case we have chosen the basis set which 
divides the given region into a finite number of non- 
overlapping elements or subregions. We make an assumption 
that whole of the density of any such subregion is concen- 
trated at the centre of that region. We can do away with 
this assumption, and not fruitlessly, at the cost of 
computational time and money. The assumption which we have 
made gives practically acceptable results with comparati- 
vely lesser (still large) computations. 

We sample our picture in the following manner : The 
picture is considered to be lying within a suitable square 
grid in such- a way that all the desired details are encom- 
passed (Pig. 3.2). The grid contains squares of equal size. 
The density of each of these squares or cell is taken to 
be concentrated at the centre. These values are taken to 
be the sample values. 
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As image processing literature confirms, when 
sampling of this type is used it is assumed that there is 
no abrupt change in the grey level of the picture. The 
levels change in a gradual manner. This allows us to 
attach one density to each cell of the grid. And this 
density, as mentioned above, is assumed to be concentrated 
at the centre of the cell. 

The above mentioned method has its implications. It 
will be much easier to illustrate the difficulty with an 
example. Consider the picture, to be discretized ,to be a 
circle of grey level one embedded in a back-ground of 
grey level 0. Following our sampling method we enclose 
this picture by a grid of size say 4x4 16 cells (Fig. 3.2). 

Fow we have to attach one density value to one cell. The 
cells which have only one density can be given that value 
but there is an ambiguity when the cell has two or more 
grey levels. In the figure cells 6,7,10,11 have their 
greater part (> 50$) filled with grey level 1 so they will 

be assigned this value. Fig. 3.3 thus gives the picture 
we are actually going to use for processing. The circle 

i 

has been converted into a square. 

It will be immature to jump to conclusions rendering 
this method useless. Let us increase the number of cells 
in the grid keeping the picture area to be same. We make 
the grid size 8x8 (64 cells). Using the above method of 
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assigning values -we see that our picture now changes to 
a cross, not at all acceptable as a circle. But when we 
make the grid denser i.e. increase the size to 16x16 
(256 cells) or more, the circle begins to take shape (Pig. 
3.5). Thus a proper grid is to be chosen depending on 
details of the picture we want to bring out. We can write 
the cell values in the form of a matrix. A nxn grid leads 
to a nxn matrix. We observe that finer the sampling grid 
the greater is the resolution. 

Thus a digital picture is defined by specifying a nxn 
matrix of grey levels. An element of this matrix is called 
a picture element or pixel. 

The unknown picture is approximated by values assigned 
to each element by the reconstruction algorithm. 

3.3 The Actual Problem 

We now have a picture digitized into STxF array of dis- 
tinct cells which are each a part of the whole picture. 

After digitization the picture is ready for taking 
projections. Obtaining projections at any angle consists in 
dividing the picture space into non -overlapping bands of 
finite, equal or varying widths at that angle. All 
the cells whose centres lie within any particular ray contri- 
bute to the raysum of that particular ray. The set of 
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values of the raysums (projection elements) then gives the 
projection at that angle (P Q ). Each projection thus has a 
finite number of non-overlapping elements b^’s. Ihe maxi- 
mum size of the projection element should be dictated by 
the presumed spatial resolution in the projection. 


We will consider projections at angles 9-j ... ©p» 
where p is the number of projections to be taken. Any 
projection gives an iT element vector e.g. b^ is the IT ele- 
ment vector corresponding to the projection at angle 9^,. 

T 

We may write b^ = (b^ -j ... b^,^) 

for each of the p projections we have such IT element vector. 
If all the p projections are considered at the same time 
then we can put all the b^'s, k = 1, ... p. , together to 
form a ITp dimensional vector. 


B = 



by defining \ tl = D( fc . 1)S4l 


... b 2(K ... b p;1 , 

k = 1 , .... p 

1 = 1, ... IT 



T 


We get 

B = (b^, ... bjj, ... b 2 jp ... “(p-lju+l "* D Wp^ 

b . corresponds to the raysum of the ray spanning region 
in S (ref. Pig. 2.1) e.g. b^ corresponds to the raysum of 
the ray spanning region CDEP in OPQR (Pig. 3.1). 
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b-j , bg ... b^ can referred to as point projections. 

These can be related to the vector x* as follows 

^(k-1)N+l = ^ x j* k = 1, ... p 

3 ' e D (k -1 )N +1 1 = 1f ** * N 

The set Ik is made up of all the cells which contribute to 

the projection b.. 

i 

Example : “b ^ 3 and I 1 51(3 corresponding set 

-D(k- 1 )N +3 is made up of all the cells which lie in the strip 
CDEE. 

To get the set D. which contributes to the projection 
element b^ we consider the following : 

Consider a cell (i, j )with center as (i-£, j--jb-) . The inten- 
sity of the cell (i,j) (cell (i- 1 )N+j)) can be written as 

x *(i- 1 )N + j‘ 

In the Eig. 3.1 09 is the main diagonal and AT is 

its projection on the axis O’z at an angle to the horizon, 
so we have 

IT = IT / 2 cos (45 - 9 k ) 

With the choice of constant ray width we divide L T into 

NIT (M > IT) equal parts. In the figure 43 = BC = CD . . . 

ST = /2 cos (45° - 0^.) and NIT = IT. If the centre of the 

cell (i,j) falls within ABOG then cell (i,j) contributes 

to b, or if it comes within BCEG then it will contribute 

k, 1 
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to the projection element b k 2 an & so on. The projection 
of the center of the cell (i,j) on the axis O’Z (at an 
angle 0^) can be written as 

^ (i~4) + (j-ik) ^ cos ^%j ~ ®k^ 

where 

q ij = tan_1 f (d-i)/(i- 4 ) 1 

Then D (k _ 1)im = { (i- 1 )F+j ; 1 < i, j < IT > 

and (i, j ) obey 

4 

( 1 - 1 ) /2 cos( 45 -O k )<[ (i-i) 2 + (j-i-) 2 l 

* cos(q i < 1 ^2 cos( 45 ° - © k ) 

After achieving familiarity with the mode of discreti- 
zing the picture and obtaining the projections we can find 
the matrix IA 1 which contains information about the ele- 
ments contributing to various projection elements. To 
illustrate the construction of matrix [A] from ITxIT picture 
matrix we take a few simple examples : „ For a 4 x 4 matrix the 


projection matrices 

for 

projections 

at 

0° 

> 

450 

, 90 

0 

are 


shown below. 
















1 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

O 

It 

O 

0 

0 

0 

0 

0 

0 

0 

1 

1 

1 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 1 

1 

1 

1 
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1 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 



0 

1 

1 

0 

0 

1 

0 

0 

1 

1 

0 

0 

1 

0 

0 

0 

A (45 

°) = 

0 

0 

0 

1 

0 

0 

1 

1 

0 

0 

1 

0 

0 

1 

1 

0 



0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

1 



1 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 



0 

1 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

A(90 

•) = 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

1 

0 



0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

1 

0 

0 

0 

1 

from 

the 

above 

matrices 

we 

can 

mak 

0 

out 

that 

there 

are 

rows 


corresponding to each of the ray of the project ion. And the 
elements ■which fall within that ray are represented by 1!s in 
the corresponding row of the matrix. The rest of the elements 
of that row are zeros. 

for example, the element b-j (first element of the pro- 
jection at 0°) gets contribution of all the first four row 
elements x-| , x^, x^ and no other elements. So the first 

four elements of the first row of A(0) are 1 ' s and the rest 
are 0 ' s.. 

With the above background and recollecting the Section 
2.4 we can define our problem more precisely as follows : 

Given : The projection vector B (= b-j, ... b^. ) corres- 
ponding to the projection angles 0-., ... © . 

1 "P 
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Derivable Information : -Is we know the angles of pro- 
jection we can find matrix A' for each 9^, i = 1, ... p 
and combine them to get the matrix A. 

Problem : With the above information and the given data 
we can form a set of equations 



The solution X, which has to be obtained using a 
projection algorithm, will give an estimate of the picture 
of whose b^'s are the projection elements. Here it has 
been assumed that we already have the analog picture and we 
are discretizing it, taking projections and reconst ructing 
the same picture again. This has been done for the purpose 
of testing our algorithms. In practice we have a number 
of projections of a real life object obtained by using 
X rays or some suitable radiations. The problem basically 
remains the same. 
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CHAPTER 4 

PROJECTION ALGORITHM (PRA) AND ORTHOGONALIZED 
PROJECTION ALGORITHM (OPRA) 

In this chapter we will deal with the iterative scheme 
we have chosen for our study. The projection algorithm (PRA) 
was found to be most suitable for our study since it allows 
incorporation of apriori data, it is applicable to space 
varying as well as space invarient problems; it is computationally 
less cumbersome and above all it is guaranteed to converge. 
Besides, like all iterative techniques PRA has a way of getting 
implemented with greater ease as compared to the other contempo- 
rary techniques. After discussing the projection algorithm, we 
present a scheme of orthogonalization of the projection equa- 
tions which contributes to its faster convergence. Binary 
orthogonalization, in particular, has been discussed in some 
detail. 

4.1 The Projection Algorithm (PRA) 

We saw in the last chapter that our problem of recon- 
struction can be reduced to a problem of solving a set of 
simultaneous linear equations iteratively in hope of getting 
better and better solution till the process converges. PRA 
does just this. 
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The projection method for solving a system of equations 
was proposed "by Kacmarz C 303 in late thirtees. Kacmarz 
prove the guaranteed convergence of this method. Tanabe [31] 
reviewed this method of solving linear equations in early 
seventies and also took up the question of convergence. 
Convergence was established by Tanabe in a complex linear 
space. Such generality is not required for our problem and we 
will assume a real linear space in our discussion. 

Once the theoretical concepts were available PEA was 
applied in the form of ART by Gordon, bender and Herman [ 5l 
Gordon and Herman [23], Gordon [21] to the problem of recon- 
struction of images from their projections. The method was 
found to be quite successful for images with varying grey 
levels. With 20 iterations of the proposed algorithm they were 

able to reconstruct the picture of 'Judy' quite satisfactorily. 

[5] . 

The system of equation that has to be solved is r 
A X = B 

where X is M dimensional column vector where M is the size 
of the picture array (N for an NxW picture). B is the 
projection vector, an dimensional column vector where p is 
the number of projections and PT is the number of projection 
element s. 
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2 

A is Wx>xW dimensional projection matrix. Thus 
we can say X e V M , eBe V lNp where and V Np stand for M 
and Bp dimensional real vector spaces of column vectors. 


An assumption will he made at this juncture which 
will hold throughout this discussion. All the elements of 
X, B and A will he assumed to he non negative numbers. 

If in the course of solving the system we get negative 
numbers they will he set to zero before proceeding further. 

M M 

A mapping f. (V *► V ) is defined which maps a M 
J 

dimensional vector into another M dimensional vector. 


f -j (x) = x 

j 


(x,a.) - b. 

— — — ol— — ■ U 

^j’-j^ 
j , ... Bp 


T T 

a. is the jth column vector of A and (y.3) = z y. 

M+Bp m 

Another mapping F (V ■+■ v ) is defined as follows 


F(i,x) = f r f 2 ... f2T (x) 
x is an M dimensional vector 


We can write this product of mappings as 

F(i,x) = f 1 (f 2 (f 5 (... fB^(x))) 

These can he applied to our assumed solution vector x 
to yield iteratively the desired solution. 
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I... we begin with an assumed solution vector x° then 
P(i,x°) will give an approximation x’ of the solution x. 
Proceeding similarly we can get a sequence of vectors 
{ x^ } , k = 0,1,2, ... defined by the recursion 

x^ + ^ = where x° is the arbitrary 

initial vector. Tanabe [311 has proved that the sequence 
x^ converges for arbitrary x° to the solution if it exists 
and is unique or to a solution x for which I |x° - x I I 
is minimum. 

The physical int erpretation of the mapping f.. and P 
is as follows: The operator f ^ projects ortho- 

gonally any vector in space onto the hyperplane defined 
by (x, ) = b... 

The operator P: projects orthogonally any 

vector successively onto the hyperplanes defined by 

(x, a . ) = b . , j = 1 , . . . $fp in the order j = Up, Np-1 , ... 
J J 

3,2,1. One complete operation by P in any order is 

called an iteration. In terms of the computer algorithm 

of PBA, one iteration updates all the elements of x(x°) 

giving x' and bring this vector closer to the actual 

solution. The operation is repeated and we get, successi- 
”12 

vely, x°, x , x , ... x n till the vector x converges to the 
solution x, if it exists, or to a vector x^ (acceptable 
solution) for which 1 1 x° - x . 1 1 is minimum. The iterations, 
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in practice, may be terminated by human interference after 
a visibly best solution has been obtained. 

The method seems to be particularly attractive because 

of its guaranteed convergence. One iteration of this algori- 

2 3 

thm requires approximately 2M * Up (= 2 it * Up = 2 I p) 

2 

multiplications and as many additions. Here F is the size 

of the picture and p is the total number of projections 

. . 2 
giving Np equations. Up << F in practical cases so the 

number of computations required, though much less than 

those required in the direct matrix inversion method, are 

still enormous for most practical images. Some way of 

reducing the total number of computations required (and 

hence computation time) in reconstructing an image is desi- 

rabl e . 

A particular scheme suggested by Ramakrishaa l 33 J in 
the context of the restoration problem, has been applied to 
the current problem which by accelerating the speed of 
convergence yields a reduction in overall reconstruction 
time. This scheme has been tested with a number of examples 
using binary and multidensity digital picture data. It has 
been confirmed that the method achieves faster convergence 
in almost all the cases though the trains vary with the 
individual pictures. 
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4.2 Aspects of Orthogonal _zat ion 

The algorithm we propose involves some degree of . 
orthogonalization and so a brief discussion of the same is 
in order. : ’ 

T 4e recall that a set X is called a linear space over a 
field K if the following conditions are satisfied: X is an 
abelian group written additively. 1 scalar multiplication 
is defined in the following manner - to every element x e X 
and each me k there is associated an element in X denoted 
by a x. 

The following properties also hold good : 
a (x+y) =a x + ay(aek, x,ye X) 

(a+fOx: = ax + Bx(a,Be k, xeX) 

( a (?) x = a (Bx) (a , Be k, xsX) 

i.e. 1.x. = x (1 is the unit element of field k) 

Consider a linear space with a system of vectors in that 

space. This space will be an IT dimensional space if it has 

atleast .W independent vectors and each of them has a unique 

IT 

representation x = y a, y.. The seminorm of any vector 

i* J J 

3=1 

in this space gives a kind of length for that vector. If 
the above system of vectors reduces to a single seminorm, 
the corresponding linear space is called normed linear space. 
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With this background we state Schmidt’s orthogona- 
lization theorem [36], The theorem states that given a 
.finite or countably finite sequence {x.} of linearly 
independent vectors of a pre Hilbert space (see def) X. 
Then we can construct an orthonormal set having the same 
cardinal numbers as the set {x^} and spanning the same 

linear space as {x). 

<J 


4*3 HllA, with Orthogonalizat ion (OPRA) 

Each of the equation in the system A X = B can be 
considered to be a hyperplane in a M dimensional space. 

If a unique solution exists then these planes intersect 
nt a common point which is the solution. Tanabe, as we 
have said before, has established convergence and the 
convergence is assured even if the equations are ordered in 
any arbitrary manner. Though convergence is independent 
of the ordering of equations the speed Is very much depen- 
dent on it and this fact has been exploited in the algo- 
rithms used. We start with an arbitrary point in the m 
dimensional space as the solution. Every projection of 
this point on a given hyperplane results in the rnov ament of 
that point towards the intersection of the hyperplanes if the 
there Is any, else towards a point which gives the minimum 
norm solution. The point is successively projected onto 
various hyperplanes (corresponding to different equations) 
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and points nearer and nearer to the solution are achieved.. 
After repeated projection of the point, it captures the 
solution point and the process is terminated. If there is 
no point of intersection then successive projection of the 
initial guess onto the hyplanes causes it to move to a point 
so that the norm of the difference hot ween the latter 

and the initial guess is the minimum. 

Motion is accelerated if the planes are wide apart, the 
greater the angle between any two planes, the faster is the 
motion towards the point of intersection or to the minimum 
norm solution. The ideal case will he when all the planes 
arc orthogonal to each other since in this case convergence 
will be achieved in a single iteration for any ordering of 
equations . 

This sort of complete orthogonalizat ion can, theoreti- 
cally, be achieved by the Gram-Schmidt orthogonalizat ion 
process 133] . However, practically such an endeavour will 
call for enormous computational effort and even after that 
wc might end up with a system completely different from the 
original one due to computational round offs. This will 
render the reconstruction completely useless as the 
reconstructed image will, in the worst case, have no 
oimilarily with the original one. 



43 


Failure or unsuit ability of complete orthogonalization 
should not be taken to mean that orthogonalization is nob 
at all useful for reconstruction procedures. ¥e can resort 
to partial orthogonalization , i.e. the more practical m-ary 
orthogonalization (m«Np). In particular, binary (m = 2) 
orthogonalization is comparatively less cumbersome to imple- 
ment and gives a saving of upto 50$ on the computational 
t ime . 

4.4 Binary Orthogonalization 

In this case instead of making each of the planes 
orthogonal to all other planes we do pairwise orthogona- 
lization. This consists in making any two successive hyper- 
planes orthogonal to each other. This is practically fea- 
sible. The following algorithm has been used in principle. 

i) Take the first equation to be fixed and orthogonalize 
the second equation w.r.t. the first equation, 

ii) Orthogonalize the third equation w.r.t. the second 
equation obt rained in (i). 

iii) Continue till the Npth equation has been orthogonal! zed 
w.r.t. the (Np - 1 )th equation which had already been 
orth ogonaliz ed w.r.t. the orthogonalized (Np-2)th 

equation 

iv) low orthogonalize the first equation w.r.t. the Npth 
equation. 
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Mathematically we may state that : 


- 3+1 = " 3+1 

" ” 3*1 


(a°, ap 



<4 sS> 





2 = 1,2, .... Up, the vectors superscripted with o 
denote the vectors which have "been modified after orthogonali- 
z at i on except in the first step when second equation is 
orthogonalizod w.r.t. the fixed first equation. 


2 3 

The number of computations required are 2MK (- 2k *Mp = 2k p) 
multiplications and an equal number of additions. This is 
same as the computations required for one iteration of PRA. 


This method has been intermingled with the projection 
algorithm in a manner such that only two equations are taken 
up into the computer core memory at any time, second is or- 
thogonally cd with respect to the first and the former is used 
for modification of the solution vector. The second equation 
is now shifted into the place of the first and the third one is 
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brought into the original place of the second and the 
operation is repeated till all the equations have been 
exhausted. This completes one iteration of the algorithm. 

I'Tot at ion ally : 

A1 

two arrays to store the two equations 
A2 

(a) denotes the contents of A array 

d = o 

Start ; 3 = 3+1 

A 1 3'th equation 

Increment : 3 = j +1 

A2 jth equation 

Orthogonalize (A2) w.r.t. (A1 ) 

(A1) - (A2) 

Use (A 1 ) for modification of x 

IF (3 = Up) STOP 

Else (Go to increment) 

This implements OPRA. 

Many of the equations in the system might be very nearly 
linearly dependent. Orthogonalization tends to seperate 
them, making thorn linearly independent and allows OPRA to 
have greater speed of convergence. 

PRA in combination with binary orthogonalization was 
implemented and the results proved our point. In most of the 
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test cases we achieved 10 - 50$saving on the computational 
time. The results are presented in the next chapter. 


■-M 

Pro jection|s 
and angles” 




Preprocessing 

block 


Reconstru- 
cted image 


Orthogonali- Projection 

zing block algorithm 



Pig 4.1 Reconstruction model 
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CHAPTER 5 

RESULTS AND C0UC1USI0W 

In this chapter we present the results obtained for the 
projection algorithm (PRA) and orthogonalized projection 
algorithm (OPRA). Various error measure have been compared, 
for the two cases, and the results evaluated. The saving 
in computational time and cost has been mentioned in each 
case. Reconstructed pictures using the line printer 
imagery have been included for visual appraisal. Curves 
have been plotted with relative errors, after multiplica- 
tion by a constant factor. The tables have entries of 
relative errors, rather than absolute errors, for a more 
effective comparison. 

1. Reconstruction of Binary Pattern 

figure 5.1(a) shows a 32x32 picture of an elephant. 
This figure was used as a test pattern for the algorithms 
discussed in Chapter 4. Results obtained for OPRA are 
being compared with those obtained for the PRA. Various 
numbers of projections were used and the results compared. 

With 10 Projections : 

figures 5.2(a) and (b) shows some reconstructions for 
10 projections. The table below gives the RMS error and 
the relative delta for a number of iterations. 
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No. of 
iterations 

RMS 

Error 

DELTA j 

I 

PRA 

OPRA 

PRA 

OPRA 

1 

5588 

4565 

2360 

1928 

2 

3993 

3432 

1686 

1449 

3 

3419 

3025 

1444 

1277 

4 

3078 

2784 

1300 

1176 

5 

2867 

2630 

1211 

1111 

! 

6 

2725 

+ 

1152 


7 

i 

! 2625 

i 

Steady 

1110 

St eady 


Table 1. 

The curves in Fig. 5.14 show "that the RMS error in the case 
of OPRA is always less than that in the case of PRA. As 
depicted in Table 1, the relative RMS error in OPRA stabili- 
zes to a particular value in 5 iterations while it takes 
PEA 7 iterations to reach that level and then only it 
settles down. The same table shows relative values of 
error measure DELTA in the two cases for several itera- 
tions. OPRA does in 5 iterations what PRA does in 7 
iterations and again the errors in OPRA keep below those m 
PRA for all iterations. This can be confirmed from the 

curves of figure 5.13. 
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Humber of mismatches between the actual picture and 
the reconstructed picture is directly related to the 
visual quality of the reconstructed picture. This measure 
tells us at how many points the two pictures differ. In 
calculating the mismatches we compared the reconstructed P 
matrix with the picture matrix 0. While all the entries of 
the C matrix were either 0 or 1 the entries in P matrices, 
after a desired number of iterations, ranged from 0.0 to 
about 1.6. All values above 0.5 were taken to be 1 and all 
below 0.5 wore taken to be 1. The values equal to 0.5 were 
also considered to be 1 . The two matrices were then compared 
to give the number of mismatches. The PRA settled down to 
an error of 3 mismatches in 6 iterations while OPRA settles 
down to 6 mismatches in 3 iterations. Thus convergence is 
achieve 1 faster in the case of OPRF though point of conver- 
gence may be different. 

The above observations with 10 projections show that 
OPRA gives a saving of about 30$ - 40$ in computational time 
and cost. 

With 4 Projections ; 

The Table 2 shows the relative root mean square error 
for four projections for the figure of an elephant. Here 
we see that the RMS error in OPRA is always less than the 
RMS error in PRA. It gives in this case an advantage of 
atleast 2 iterations. See Figure 5.2(c) and (d). 
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TFe 

Number of 
iteration 

MS 

DELTA 

PRA 

OPRA 

PRA 

OPRA 

1 

6190 

5793 

2615 

2447 

2 

5340 

4803 

2253 

2029 

3 

4795 

4412 

2025 

1864 

4 

4439 

4176 

1875 

1764 

— 1 
5 

4377 

4017 

1849 

1697 

6 

4314 

- 

1822 

1648 

7 

4052 

- 

1712 

- 

8 

L 

- 

- 

1623 

j 


Table 2. 

The same table show relative DELTA. This again veri- 
fies that error is lower in the case of OPRA. PRA reaches 
to a comparable value two iterations after the OPRA has 
reached it. The convergence is thus faster in OPRA and we 
have a saving of around 305* in the computational time and 
cost . 

Mismatches were calculated for 4 projections and a 
number of iterations. PRA was seen to settle down in 7 
iterations while OPRA settles to nearly the same reconstru- 
ction in 5 iterations. We again have a saving of about 30$. 
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2. Mult ideas it y Pictures 

a) First we consider the test picture shown in Fig. 5.1 (b) 
Reconstruction with 10 projections : 

To aid visual comparison mismatches were evaluated in 
the following manner. All densities, in the reconstructed 
picture, above and equal to 0.5 were taken to he 1.0 and 
all those below 0.5 were taken to be 0.0. Mismatches were cal- 
culated for PRA and OPRA and the results are given in the 
Table 3 below. 


No. of 

Mismatches 

RMS 

iterations 

PRA 

OPRA 

PRA 

OPRA 

1 

99 

97 

6694 

5628 

2 

40 

40 

5059 

4304 

3 

33 

32 

4509 

3924 

4 

29 

27 

4275 

3751 

5 

26 

25 

4136 

3646 

6 

24 

22 i 

\ 

4034 

3570 

7 

22 

17 

3951 

3509 

8 

20 

17 

3883 

3458 

9 

18 

17 

3824 

3415 

. 3 ! 

17 

17 

3773 

3378 

— 


Table 3. 


I.i. ». 
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We see that both the methods start off with about the 
same number of mismatches but OPRA converges to 17 mis- 
matches in 7 iterations while PRA keeps on converging 
slowly and reaches that solution in 10 iterations. Pig. 
5.3(c), (d) are included for reference. A clear saving of 
30$ is observed. 

The relative mean relative error in the two cases is 
plotted the Pig, 5.15 which shows that OPRA, has lower error 
in general. . The same holds for the relative RMS in two 
cases. See Table 3 and Pig. 5.16. 

Reconstruction with 8 projections : 

Curves in Pigure 5.17, shows a plot of mismatches 
versus iteration. Both PRA and OPRA have same number of 
mismatches in the initial steps but as more iterations 
are performed OPRA accelerates to the solution. OPRA reaches 
27 mismatches in 9 iterations while PRA reaches the same 
value in 13 iterations. The it erationwise mismatches 
are shown in Table 4. Again a saving of about 30$ in compu- 
tational time is thus achieved. 


Mis- 

No. of 

1 

2 

3 4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

matches 

it ora- 
tions 














PRA 

118 

63 

51 41 

48 

35 

33 

32 

28 

28 

29 

28 

27 


OPRA 

118 

61 

52 38 

36 

35 

31 

30 

27 

27 

27 

27 

27 


Table 4 




Table 5. 

The relative MS errors are showi in the Table 5. The 
observation confirms the results obtained above. 


Reconstruction with 6 projections : (Pig. 5.5(c) and (d)) 


No. of 
itera- 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

tions 

PRA 130 

89 

82 

78 

76 

74 

70 

65 

65 

63 

60 

59 

Mis- 

matches 

OPRA 128 

83 

77 

75 

71 

69 

68 

64 

61 

59 

59 

59 



Table 7 
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The mismatches, in Table 6, show a gain of about 20$ 
and the table of relative RMS error (Table 7) shows that the 
error in the case of OPRA is always lower* see Pig. 5.18. 

Reconstruction with 4 projections : 

With 4 projections we get in 3 iterations of OPRA what 
wo got in 5 iterations of OPRA. A saving of 40$ is achieved 
in this case. See figures 5.6(c) and (d). 

Prom the above discussion we see that for this test 
pattern, for all sets of projections 4,6,8,10 we achieve 
faster convergence with OPRA and there is a saving of 20$ - 
40$ in computational time and cost. 


b) A three density pattern 

In Pig. 5.1(c) we show a pattern having large trian- 
gular areas of densities 0.53 and 0.93 and small areas of 
density 0,0. Reconstruction for various projections are 
described below. 

Reconstruction with 6 projections : (Pig. 5.7) 

Mismatches were calculated by fixing some range for the 
various densities after observing the P matrix for several 
iterations. Values of density below 0.15 were taken to be 0. 
Values between 0.15 and 0.67, including the lower limit, were 
taken to be 0.53 and values above 0.67 were taken to be 0,93. 
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The mismatches obtained for 6 projections in the two cases 
are shown below. 


Wo. of 
it era- 


1 

2 

3 

4 

5 

6 

7 

8 9 

10 

11 

12 

13 

..ti.oaa.-_ 














Mis- 

PRA 

299 

155 

116 

95 

81 

75 

67 

62 58 

54 

53 

54 

53 

OPRA 

I 

285 

i 

157 

114 

96 

84 

76 

71 

67 60 

57 

51 

51 

51 

matches 

i 


Table 8. 

The mismatches show a behaviour different from those in the 
test picture 1 (elephant) and test picture 2 (windows). Ini- 
tially PRA performs better but after a few iterations OPRA 
catches up and converges to the solution much before PRA. OPRA 
converges to a solution having approximately 51 mismatches in 
11 iterations while PRA tehes 14 iterations to settle down to a 
steady solution. Gain in this case is nearly 21$. 

later the picture was reconstructed using a number of 
projections 4,5,6, 9 to study mean relative error for com- 
porison. The error was plotted for all these projections with 
a fixed number of iterations- (3). The curves are shown in 
Pig. 5.19. Wo can see that the error remains low in OPRA. We 
can also see that increasing the number of projections above 
6 is more fruitful in the case of OPRA. 
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c ) G-r a d ed s quare 

The picture shown in Pig. 5.1(d) consists of four 

slices of different densities, namely, 0,0, 0 . 42 , 0.56 and 

1 . 0 . 

The test picture was reconstructed using a number of 
projections. Representative results are given helow. 

Reconstruction with 5 projections : 

The test picture was reconstructed using 5 projections 
(Pig. 5.11) and varying number of iterations. Various errors 
wore* calculated. The way the mean relative error varies with 
the iterations is shown in the Pig. 5.20. At times OPRA 
performs almost same as PRA or slightly worse but ultimately 
OPRA gains over PRA and settles down to a t steady error 
after 11 iterations while PRA settles down after 14 or 15 
iterations. A gain of over 20$ in computational time and 
cost is thus obtained. 


RMS error in the two cases have the following trend : 


No. of 
it ora- 
tions 


3 

5 

8 

10 

12 


PRA 

7316 

3115 

2866 

2737 

2693 

2664 


7057 

3085 

2866 

2737 

26,88 

2656 

SOI 


Table 9 
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OPRA starts on a better note but from 5 to 8 iterations its 
performance is no bettor than PRA. After 8 iterations 
again improves over PRA, giving lesser SMS error. 


The overall nearness factor, DEIff A, has a similar 
Reconstruction with 8 projections : 


No. of 
it ora- 
tions 

RMS 

error 

Mean relative 
error 

DELTA 

PRA 

OPRA 

PRA 

OPRA 

PRA 

OPRA 

1 

6264 

5763 

2250 

2241 

1679 

1603 

2 

3894 

3558 

1370 

1351 

1044 

990 

3 

2886 

2638 

978 

962 

774 

734 ' 

4 

2397 

2179 

805 

789 

642 

606 

5 

2105 

1908 1 

701 

685 

564 

531 

6 

1906 

1726 

632 

618 

511 

480 

7 

1766 

1596 

586 

571 

473 

444 

8 

1662 

1500 

552 

538 

445 

417 

9 

1582 

1427 

526 

513 

424 

397 

10 

1519 



1369 

J 


493 

i 

j 

407 

381 

1 


Table 10. 

The above table (and the curve in Pig. 5.21) show the relative 
RMS error, mean relative error and DELTA in the two cases. 
Once again in the case of OPRA the error values are always 
lower. Wc gain here atleast two iterations in the steady 


state 
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Reconstruction with 11 projections : 

With 11 projections the following error values were obtained. 


No. of 
it era- 

RMS 

error 

Mean relative 
error 

DELTA 

tions 

PRA 

OPRA 

PRA 

OPRA 

PRA 

OPRA 

1 

7217 

6484 

2329 

2327 

1934 

1804 

2 

4491 

3466 

• 15H 

1289 

1024 . 

959 

3 

3735 

2836 

1266 

1063 

1001 

789 

4 

3443 

2549 

1160 

951 

923 

709 

5 

3303 

2387 

1106 

886 

885 

664 

6 

3217 

2282 

1072 

844 

862 

635 

7 

3157 

2206 

1046 

814 

846 

- 

8 

3113 

2147 

1026 

790 

- 

- 

9 

3078 

2098 

1010 

771 

- 

- 

10 

3050 

2060 

997 

- 

\ 



Table 11. 


All the three error criteria have been evaluated and compared. 
OPRA performs except ionally well in this case. AL1 the errors, 

RMS, mean relative error and DEITA keep much below those in PRA. 
There is a saving of over 50$. The curves can be seen in Pig. 

5 . 22 . 

Por this test picture we see that with 8 projections we 
are performing better as compared to 11 projections. This can 
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happen and is a normal phenomenon 26 . This can he explained 
by observing that more the number of points on the edge of the 
ray , the greater is the error since the point can be assigned 
to only one ray while the square, of whose this point is the 
centroid, might be contributing almost equally to two rays. 

The gains with OPRA have so far been stated in the form of 
$ gains. Wo can get the idea of savings in terms of the units 
of time and money by the following : 10 minutes of CPU time cost 
around Rs. 120 on DEC 1090. Thus if PRA takes 10 minutes to 
reconstruct a picture matrix then a saving of 30$ on the computa- 
tional time will mean a saving of about 3 minutes and in cost of 
about Rs. 36. 

Conclusion ; 

Eror the discussion above and the conclusions drawn at each 
step we can say that binary orthogonalization does speed up 
the algorithm (PRA) under study. The four test pictures that we 
have taken for our study, provide enough proof of the superiority 

of OPRA over PRA. 
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APPMDIX k 

A. 1 Reduction of Dimensionality of the Problem 

A three dimensional object can be considered to be 
made up of several 2-D slices put together. If, for 
example, f(x,y,z) is the 3-D object then f(a,y,z) represents 
a slice at x = a. These 2-D slices can then be reconstru- 
cted using the projection algorithms proposed in the study 
and they can then be stacked up to reproduce the original 
3-D object. Similarly an R dimensional object f (x-j jXg, . . .x^) 
can be reduced to slices of (N-1 ) dimensions, viz. 

:f(n ,Xp,Xj, ... x n ). The (N-1 ) dimension problem can be 
reduced to (H— 2 ) dimension problem and finally, using such 
reductions -we can get N dimensional problem reduced to a 
number of 2-D problem. 

We favour reduction of dimensionality since it redu- 
ces storage requirement and simplifies the algorithms. 

A. 2 Reduction of Storage Requirements 

In the reconstruction problem we have to deal with 
large matrices. Even a 4x4 grid produces an A matrix of the 
order ,Npx16 where Np is the number of projects. For only 
4 projections the size becomes 16x161 But we observe that 
the matrix A has either 0's or I’s and most of the elements 
are Os. Such matrices are called sparse matrices. In such 
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circumstances storing up of whole of the matrix will "be 
unwise. There are a number of standard procedures which 
can he used for storage and retrieval of such matrices. 

In one procedure we store the row numbers and column numbers 
ol the 1'S3. This reduces the storage requirement (as in 
1(45) Oh. 3) upto 50$. Also we observe that there is a 
sort of circular symmetry in that the last row can be 
obtained by right shifting the first row a fixed number of 
times, the last but one row can be obtained by similarly 
shifting the second row and so on. This means that it is 
required to store only half the number of rows of A. 

This presumes that the number of projection elements, in 
any projection, is even. A saving of 50$ is again obtained. 

A. 3 Compactness of Figures in Digital Pictures 

The square of the perimeter of a figure divided by its 

2 

area, p /A, is a classical measure of non compactness of 

the figure. In plane geometry p is the normal perimeter 

2 

of the figure and p /A takes a value greater than or equal 
to 4tf, 4i* being the lower limit which a figure can take 
only If its a circle, .ill the other figures are less 
compact than circles. For example, for a figure of a square 
on side s , * p « 45, A = 5, P 2 /A = 16 > 4*. But for digital 
pictures certain squares or octagons can yield smaller 
values of p 2 /A than do the digitized circles. This depends 
on how the p Is measured. But there is no family of shapes 
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2 

lor all ol which p /A takes on a unique minimum "value • 

Digital Perimeter 

There are two different ways of measuring p of a 
figure s which require attention. 

1) P can he taken to be the area of the border of 
where the border of s consists of those points of s that 
have horizontal or vertical neighbours not in s. 

2) P can be evaluated by applying the border following 

algorithm to s [51] [52] by counting 1 for each 

horizontal and vertical move and /2 for each diagonal move. 

The method 1 will be taken to give perimeter P-j and 
the method 2 will be taken to give perimeter P 2 . 

Exl 

1 1 1 

1 1 1 

1 1 1 

a digital square ; for this P 1 = 8 , P 2 =8 (since all 
the 8 moves are in horizontal or vertical directions) 

Ex. 2 Ex. 2. 1 

1 1 1 
11111 
1 1 1 
1 

a digital diamond : for this P-| = 8, P 2 = 8/2 (since in the 
case of P 2 all the 8 moves are diagonal) 
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Ex. 1 

1 

1 

1 

a figure of width I s for this P 1 = 4, P 2 = 6 (since 
in the "border algorithm we start and end at the same point) 

2 

Digital compactness ratios, P /A 

Consider a regular octagon 0(h,k) whose horizontal 
and vertical sides have h points each and the diagonal sides 

have d points each. 

Then P 1 = 4(h-1) + 4(d-1) 

P 2 = 4 (h-1 ) + 4 (d— 1 ) / 2 
A = (h+2 (d-1)) 2 + 4 |(d) (d-1) 


h + 4h(d-1 ) + 2(d-1) (d-2) 


Wo observe : 


1) Dor a square d = 1 


So. 


p. ® 4(h-1) + 4 (d-1 ) = 4 (h-1 ) 


P 2 a 4 (h-1 ) + 4 (d-1 ) /2 = 4(h-1) 


ii sat h 

* 

_ 2 , . „2,. 1 6 (h— 1 )‘ 

Thus P^/A sa P 2 /A « ^2 — 


16 for large h 



64 

2) For a diamond, h = 1 
There, F 1 = 4(d-1 ) 

P 2 = 4 (d— 1 ) /2 
= 2d(d-1) + 1 

Thus 2 ip A = P 2 /A = 32(4-1 ) 2 /[2d(d-1) + 1) 
for lnrgo d ip A - 8 and P 2 /A - 16. 

3) For an octagon with h = d 

P 1 = 4 (h— 1 ) + 4 (d-1 ) = 8(h-1 ) 

P 2 * 4(h-1 ) + 4(d-1 ) /2 = 4 (h-1 ) (1+ 2) 

A = 7 h 2 - lOh + 4 

ip A . 64(h-1) 2 /(7h 2 - lOh + 4) 

2 

tor largo h P-/A “ 64/7 = 9.14 

P 2 /A = 16(3+2/2) (h-1 ) 2 /(7h 2 -l0h+4) 
for large h P 2 /A a 16(3+2 2)/7 = 13.32 

4) A digital circle of radius r can he defined as a set 
of points whose distances from the centre are Hess than 
(r+*), 

2 

P^/A for such a circle was found to he nearly 10.2. 

p 

P 2 /A was about 13.7. 

2 2 
Thus P j/A is smallest for a diamond and P 2 /A is smallest for 

an octagon. 
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jiPPBNDIX B 

B.1 Numb or of Views Required for Reconstruction 

If a reconstructed image is to be uniformly resolved 
to a resolution d of a completely, unsymmetrical object, the 
numbers of discrete views must be at least HU H2J 

n = irD/d, where D is the dimension of the object 

Thus for a resolution of 2 cms. in imaging a head 20 cm 

in diameter, we need nearly 32 views. 

But here we have assumed that the object is completely 
un symmetrical and random but this is not the most represen- 
tative case. The practical objects of interest, be it any 
part or. organ of the human body or be it a celetial. body , 
there is a great departure from complete randomness. This 
allows 5 to 20 projections to suffice in usual reconstru- 
ctions. 

B.2 Soperation between Projections - limits on the Single 

of Projection 

If we take projection at 0° (shown by broken lines) 
then for the first projection element OB the pixels 

involved will be x^ , x 3 > x 4 * 
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of projection 

If wo ino to no o the angle slightly then we still have only 
, X£» Xji x^ contributing to the first projection element. 
Wow v/o go on increasing the angle till the first ray just 
min nos the pixel x^. At this instant the angle has been 
named Wow these two projections (at 0° and at 0 -qjj,) 

will have some different information to convey. The two 
equations corresponding to the two projection elements 
OB' and AD will not be totally linearly dependent. 

The situation depicted in the figure B. 2,1. refers to 
a difference of one cell width (07 - R5) between the 
intersected lengths at the point of entering s and the 
intersected length at the point to leaving (7) the object. 

w© have to find the angle ^y which any two 
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projections must be s operated to cariy different infor- 
mat ions . 


From the right angled triangle VOT 

ov 

tan Q bif = m (b. 2.1) 

From the right angled triangle SRI 

tan © b -jji » jjjr (B. 2.2) 


From triangle BAG 


AB » /?., coo (45 - 0j)j F ) since AT = /2 N cos(45°-©- D -j- p ) 

no AB /2 cos (45° - Qp ip ) 

cos' SpxF ** c °s 0 BIP 


« 1 + t an Ojj-j-ji 

07 *= BC « 1 + tan 
From (B. 2.1) and (B. 2.3) 


OT a 

From (B. 2.2) 


07 

tan © DIp 
and (B* 2.3) 


+ 1/tan ©j)j P 


RS » RT tan Q mF 

m (OT - OR) tan ©pjp 
« OT tan Q dip - R tan ©pjp 
« 1 + tan 0 BIp - $ tan © DIP 


(B.2.3) 
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Since OV - RS = 1 

1 + tan Qjj-j-p - 1 - iaii ^DIP + ^ ^ an ®DIE = 1 
N *fc sn ®j) jp ” ^ 

0 D Ip = tan" 1 (1/N) 

for N = 4, e rip = 13.5 so the projections must he spaced 
atleast 13.5° apart to give useful information. 

B.3 Error Measures 

The following error measures have "been used to compare the 
algorithms under test. 

1. The overall nearness factor DELTA, 

DELTA = ((1/N 2 ) I (0(1) - E(I) 2 )^ ; I = 1 , . . . F 

I 

2. The Root mean square (RMS) error, 

RMS = (I (0(1) - Pd) 2 /^ (0(1) ~ AVG) 2 )^ 

I = 1, ... N 2 AVG = (l C(I))/R 2 
5 I 

3. The mean relative error (MRERR) 

MRERR = ( E I 0(1) - P(I) I ) ^ 0(1) ? 1 = 4, . . . F 
In all the above definitions 0(1) are the picture matrix elements 
and F(I) are the reconstructed matrix elements. 
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APPENDIX C 


The hasic code employed, with the approximate density values 
on a normalized scale, are as follows : 


0.00 


0.15 

= 0.22 

+ 0.25 

) 0.29 

1 0.33 

Z 0.37 

X 0.40 

A 0.42 

M 0.45 

0 - 0.53 

0 = 0.56 

0 + 0.60 

0 •+ , 0.64 

0 + , . 0.67 

0 + , . = 0.79 

0I,i- 0 • 85 

OX, .HO 0.89 

OX, . H B 0.93 

0 X , . H B V 0.97 

0 X , . H B V A 1.00 
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Pifr.5.2 t Reconstruction of the picture of the elephant s U?Ra.lO,ilEQe473 
JIT.6, (a) PRA, (b) am* ’ OT«3| HPR«4, NECtl 89,3*17*7, 

( 0 ) PRA, (d) OPRAt tfPR.8,ii£§ vJntQ»378,HIT«5, (e)PRA,(f) OPRA, 
NPR»8,i^Q«378,OT«5. 
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* Pifl*^*** f"' 1 & 

* • naeaftiM 7 v ^' r 

fc «i Vi t v 

,‘“g > i i i ^ 

* **«a h ^7 

t f ) -‘MU’/ « 77 V 

* - ♦S))r( ,, S 4 ? 7 V\ M 

* --if 1 Y'f j v y a ’• 


-- 1 *« 4 1 7 - 




f 
* 
£ 
* 

■ *• 

fc 
* 
* 

:*» * 

V x *' '< ■* M Y|)| 9 ^ * 

-*V 99 $* 

4 ) * 1 1 1 ? 7 fUtH 4 

mm r * 

Y > M ?.? V'Ytttf ♦ 
+ *111 7 v Afttttf * 

Vi 1 7777 ^^^^ 

MU 7 r<Atp#»* 
HM’? v^*| 9 « 


f - 
l r- 
1 ♦ 25 
1 1 1 

7 7 

}* mm m m 

1 ****.«**, ~ <*«. *« 


wia#» 

A M t 

y v ■* 

v ' 

* : 


* i * * * * * 4 t : + * * * <f t * * * 4 4 * * 4 ' * + * * * & 4 & 4 * t 

(b) 


*********** **V **** ***** * ** * *** 


********* t*Uittt***n*i»*(i*f Hfri! 


♦iiiiiofffffifWMZ)— * 

♦ifii#8»»i*#ff*AXZ)- 
♦niiifHfffffiiMMAZis * 

♦» 90 wiiitifffAZ))nii— * 

♦if ■ •■••••••o ) ♦ > ) j i zx»— - * 

♦f®®*l®iii®MA3aa+)lZXa-x— a. a 

♦ XxTz8XZ)U0MM)nzjcMM0)))J) + + *MAZ + * 

♦ x 22225 f? 2 MM *®* zzzxx * zz #»» w * 
t ,«VY* 2 A?’j«? 1 ? 5 fJ M 2 alZZZZXXX1 *» ,# * 

MMMMMMMH ) ) ) j ZXMMt ) ) ) ) 1 1 i •••> < 
- i XZ1 XAHH1 1ZZXAAM4 ) ) ) } r/.ilil * 
) ♦ t ZM8ZZXXHHMH4) ) l UZZiffS* 
-»*4lZMl AMMMHftOl) UZZXXiffS* 
— + )X4ZM8if«*2ZZlXXAAff»»» 

SS at w 

♦ 4 4 *****«4*4*«****** 

(c) 


♦ aafiyftiifva 1 ' t- - iy) * 

♦iiiiiffiffiiAVAZzn t 4 

♦ i#ii|i«9Wii«y?M n 1 1 — * 4 

♦ •■•if 9 *i»i»») ))nr/xa — - * 

♦ 4 ’ll \ 77 — S-— S. + 


♦ ft* t7,«*7) ^nnnt 1 1 7r ‘«mm> 

♦ 9 *»B»iXXA»^aiiZX 7 ft v'.' V 7 « 4 ia. * 

♦ - ‘ OBB'' *< ) 7.?. A r ' «»87ZT7 7.y 

♦ AMl* 017 f'»N )M||I! 8 «ac* 

1 M»rn> t $j 8 i)}ft* 

♦ ) 1 1 ?v *.••’) m i < •/'»;»*•»* 

♦ «t »:■? ! XMCI4M 1 7,7.y AA ’! + ) hh ' ? x*i»a* 

♦ :=ni Z"f*77.v 1 17.7.V »9|)5S4 

♦ --= 4-1 v 'M y a ■« i ) 1 1 7 Y”n*@f|v* 

♦ — —h a* r«M»a«*n 7 i ?y. a A'sai’’ * 




Pi‘% 5.6 


Hoooantructioa of vindows : 1IPR»4 ,^0=199,^11^2, {a.) PRA,(b) OPRA; 

t jpR« 4 ,JriSQai 89 , (c)pra,,f^^ 5 » (a) opra,^^. 


********************************** 


********************************** 


*#OOMMZi ♦■■■♦Izzt • ) | XZXrtMOO* 

*08808 AZ1 ))OZZZ— -a**} UZXAXMO80* 
*8MMMM8 MMMZXXj A •••♦) J tZXAMMMMMQ8* 
♦WMAMOMOSOOOWM ••«>iZXXMMMMOQO§M* 
*X A MO8O00080O-B--BB 1 2X AMMQ0O*UM»»n * 

samamin s^fe' 

«g 8«8inn: a8iiiiaiiMg; 

gS«Sis8K^~:™Hi!SSSg: 

♦8808XMXX— ••*••••••— — «♦) ZXMM9* 

♦OMUM)- • • ■••MSaatSS) 1 ZAMM* 

*9ZX" - — a*a*s«aa+) IX AM* 

***** <:***************«************ 

Ca) 



♦80M88M8AAZZX1 X — ♦)) 1 ZXAMMWMM08* 
*HMMf»0M0n000H« ••«) tZXXMMM»0008»* 
*AAM880»iM*8»0aa-«aalZAAM»8088»0®0* 
*ZXAM8f9V«M«t *f »*) 1ZMM8989®0»0©»«* 



*MMM«00880989®MA0MA88«0XXAAAMMMMM* 
*A888a88B888888AMMAXAs*lZZMMMrtM8H* 
*M880888888888MM8 b**«3c)ZXAMHH 88B* 
♦M8888898B8MMH I ) ••••—a* 1 ZXMMM88* 
*8flM#«8®MMYA*»— *— — — a*)Z»MM88* 
♦8M88A8MM— —————a* )ZXMM»* 
♦OOMOZt* - • a.»*«aa—*«a) JZXMM* 

*880«-» •••*♦*««»**♦ ) t XAM* 

********************************** 

(b) 




********************************** 


*9®8ZX) 

*a©f©MA„ 

*8«f888AZ**aa 
*ifl89998MZl*a 




®MZl*ah asasl) 1ZA899999 
899MZt>* -aalXAMB 
898AZ <-*X899i 

8f*tr^U98itif999 

i 9 mSSU 2 |Ui 

fsspt:?;2nifnf 


9088MXU) *♦«+ )•♦«•«« ) X8 



SHMasgiinssre 

a 

********************************** 

(O) 


********************************** 

*89MZX)J)->+)ZlaZ)— ^•♦♦J)))X8Z88* 
*8089MAZ+4aa*l)lZa-a**!ZZZXHM8888* 
*898880AZ)*++lX)+a«**itZZAM8888§p* 



SSsisssI 
ttpsHsnis 

^■■■Si099BMAX**a *al lMtsiofilil 
^■9909 88MX 1 1 ) a*a-a-a* 1 A099ilfil 


M :sggfi3s::s|»usuus: 

00998H9899099009O00990* 

f 

* 
* 
* 
Jft 

MA1Z) ♦♦»♦)«♦*•*») ZM00§f M * 
•098 AZ1 1 )♦) )a l ♦♦««* ) 7898809* 
_990AZ ))♦♦♦) at t ) }a*a*ZA80®0* 

*0M8Hssas))* + )as)alZn) + )a)1XM80* 

********************************** 


(d) 




■ . 


Pin. 5. 7 l Reconstruction of the pattern} iJPRs:6,^iiQ=24 ? II'P=1 ,(a.)PRA, 
(bj OPM; iTPR*6 , JOQ=2S4 , ( c ) PRA, NIT=14, <d) OPRA, 
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**{*>::’ A* 7. t i l' , <YX'V ,, !4»>'«* 

fUyon^_M?t > »? )1-S4ar».XAAA ^’«yn iM'.JO-* 


f 9 »® 8 **AXZf HZ 7*11 ixx‘ft*i 


■«***•* 


*9#«9a9*A AXAK '57 1 y x I y. **?..?* 

'<■*?.* x f x n 

*,g£«f 7 . n x 1 x x a n x .v* 9 »»t ®»99 * 
*i§® 99 ® 9 »® 89 *w 8 > *?a f a ■* a A ••*»a* 9 Pf ® 9 ® * 

*999999099909 «»*?•<••■, t - ^ 

*8»9®999£9Pi99 M 7M7X* - «!^*®«9*9tl9»* 
*9®9®99«9f®9t*«/,^iri»»0*|i*****®g5ll»* 
* 9 f 99 ®»* 99 P» 8 ft* nz) x ’’apatt »•«••#« v * 

f »©•»»• »io*< xy i )s))8i»n «^piD*«Sar* 

t§f 99i9®8i« a x 7 1 )a+ts 1 1 7Aj’.‘*a9>a»>#aa** 

*if»aa«M»x7t t ))s+*a)) 

*§9999*^X1 1 Jl > *- + + «»* 1 17iA"i’«H00»09* 
♦90HAAXZ! )U*4s»as«ni 1 v.Y'-.M^paju^d* 
*977,7.1. 1 )+tS8!SS««ss«+ 1 t 1 7,Y V’^aHttHta* 

******* +* i* 4 ** + ** ***************** 

(a) 


********************************** 
♦HfMMAAZl ))♦♦«♦• **8>tZZXXAMM9989* 

!2222!2^ J< 2n? ), l* s+slzxXAAMM ^ noo o f) * 

*iftt»#MXZ71 ZXZa ) t ) X A AM«MMM8M99§»* 
♦••i»9B«MAXMMM7l XXI XAMWMM«MM«9§if* 
*9M9999B«00ftMMZAAIAAA«*MM88B9«5** 
**t|t*f»»»*M8§XAAl AAMAX *88909909* 


♦tf*|9»»gf«i®a?<zsA7.AM«e»*a*#®ff®i 

*tiM«»Bfi 9 »«»axMMABa» 99 « 9 «B«iiii 

*i«»*»aatat e«A 7 ) zzt xMaaaaaaataiii 
♦•••mwmmxzi * 1 1 )ZXAaoaa»a»»» 9 fl 
♦••••••••wmaxzi s) } * 1 7.AAMaa«aaa»»f 

*fia»$#e9MAxzns*5anr — 


iZAMMaaoaaaai 


' — tt # ■— * ♦ «n rr 1 wpawaaaw -*■ 

*999999 m m a z i i n a ♦ + * ) t zz a mm»9»8999 * 
*• aaa » m a x z u ) ) ) - ♦ + * ) > 1 7 , x a Mawaaaaa * 

*i»o«Axzi i )+**=-s+sj)i7zAMMw«aa*a* 
*aMMXZl ) ) +s!!sa»it»t } i 1 ZXAMM«wa®a* 

********************************** 

00 




********************************** 

»aaMxi!*+zt Azz , »*f*»*)zi**4)zzzi*»f* 
Y*amaxzi)*)zAi )++*-f)z«*)ZAMMeiii* 
*e3S;R?I9A) >44111 *4 4 - ) 4 se*ZXA 0 «® 998 * 
* 9999999 x 1 * t*az* )♦)«-» 1 t am&«9«999* 
19999*814 1 4 4 - ) ^ ) XZ 4 4 1 Z 7 .MB 0 « 9 «ir ‘ ‘ 


*01 

♦Bf 
Jf 

;bl 

» 

♦iMtaBrtMMAAxzz+zzuxjzxiMftaioiiaY 
♦ifaaaMHMzzx l j t ♦ ) ) *mx— »« i naaaaa* 

♦S®»M.MXA «)1 Z— — + > ) 1 1 )* 4 ) 1 A 999 * 

♦aaxzzt— **-—*-***) t XMXMA 7 .) t xaai* 

3 3 

********************************** 

(<0 



♦ximaBaaa 

•aMmmaaaaaai 

~ 0 § 9 X- 1 ZXttBVI 

( 91 *A AZatai 

(98) 1119991 

iitMBaxMtSff 

iftfAMMizMaifi 
, 9998 AMA 1 ZZZM 9 
99 MM«MiAA)ZZlZH 9 



********************************** 

♦a®K X 1 1 S 4£5 A 7 /,) 3 + 3 S 4 ) Z * a ♦ + 1 Z 7 7.999 * 
fa«®BXZl )♦ >ZA 1 ) + + 3 »+)Zs*l 7 .AMM 99 a 9 * 
4999990 A 1 ) ♦till 3-M--)+sc+ZXA»«*00f * 
*99 99999X 1 4+ssZs ) * ) 3 - 3 1 1 AMP000099 * 


Sli 


| 9 | 909 iMa 9 XMa 999 j 

lM99999A«MlZM*aa| 

*ii 9 | 83 ® 9 «®» 989 A MAI ZZZMSiiL 

4 |iii 9 « 9898 MMeM 1 AAUZI 7 .M 00 f 0 »*#i£ * 

*1®9900 H M m A A X Z Z + 7 Z ) A X 1 7 X 1 M8a0«9®0 * 
♦iittStMMMZZXI ) 1 4 ) 5 ♦MX— C31M909984 
*ii®»9MAAXl )1 Z1 »* -344S-34 1M9999* 
**9®MMXAA)1Z—- — - 4> ) 1 ns4)1A*»®4 
*99X7.7,4— **—3-ss*)t XMxVAZ) 1 AHB9* 
* * 
********* ************************* 
(a) 


Pin:. 



5.8 1 Reconstruction of the pattern t ICP1M, tfEQ=l89,NT**3, (a) PRA, 
(b) CPRA; ilPit= 5 , »I 3 Q s .S 32 f SIT« 3 f (c) PRA, (a) OPRA. 


********** ************************ 

♦08AZX) ))»♦♦) 11 4Z1 «—«)))) mxx#«* 
♦980MAX 1 4 44 44l ) ) x- **4 ) ZZXXAMM000* 

1 00998MX1 444 41 A 1 «— «4 ) 1 ZXM*»f»*«H»i* 
•••••80 A 1 1 ) ) ZZ*~« ) ) ) ZMM0808888** 
*888888888 A AX 1 4 -*«1 XAM0909999909* 

l 888 |S 8 «SSM|»SIS| 8 SSSS«Se * 

I888MAA41- ••*41888888888* 
__ mm J 88AM)«— •—♦188888888* 
|iiBi8§flMZX*-*"**»*« - "»4 1 A0}i)0§* 
iBlfltMMl |***«**4* ) 4«*s»4lM88iii * 
}9B9AZ 44««4 ) 44-44* l ) ) 4*44)Z«88*§* 

iM84*4S"»*a*s+)s+)sizi 1)4)4128888* 

3 S 

********************************** 

(a) 


********************************** 

♦08MZX) ) )«4*)Z1 47 !■•**) ) I ) 1 17-AX08* 
♦0080AM 444441 ))X- *a4)Z7XXMHM980* 
♦89888MX ! ) 444 1 A 1 s»o*4 1 ] ZA M88888884 

*9889890* 1 U ) 2 1 *— ) )I XM808888888* 
♦8899 888908 AX 14 *a*l X AM#8»00i*»00* 

!s:s!:::s»s»»:ss»»siu::»s|i 

♦■999900998998989 9M09999080998969 ♦ 

:s!i8::sg;5i5:355T?s?5;s ::;:;::s: : 

♦998989990MA2Z) *«• -••■4 1 M***ti§f* 



♦0000M A )♦**♦) 4 4* *4*1 1 ) )*♦♦ )Z9fJ000* 
*989) )••••*♦*♦) ■4)*2Z11 i)l4lZM888* 
* ss 


********************************** 

(b) 


• • V 


** * 8 ¥ *** #♦***.*#** ***♦♦***# ******** 

♦ *9 X.7. A Z0 A 07. ■ 04 ) ) 47.X 1 X A X7 1 7.Z7 1 7X808* 
*998907. X » M7, ) t ) 1 * ) 2) X 1 Z7.X7, ) 11 H8888 ♦ 
♦908090XXXK)71)«1 11 AAAAZ4)7.008808* 
♦8889009M Z 1 X A Z7.« t Z 1 A AM 1 X X A99808 88* 
*808098 9881 22ZA41Z>ZXZ7.H«0808 ••••* 

♦ • V09899990MX ZZ l A A? A) Z90900998999 * 

♦SiSiSiivSSssifizSfMoisssSiSiiiSvt* 

♦•9w30000«i08i9ZMX*a099i0i0«99»tt9* 

♦•Iff - 

♦ 0981 


**®0ii0000fi000M8M*.X v 0i0i980f0890* 
♦•*99099909089* 1 M A A A t * 1 *8*99® 0981 
'•9888***«90AA7,)7« 44*4)0990f 
k008»***W«HMMA) -I—*) 4 1 XXH88I 


I* 

♦0»8f****0M8 A A7, ) 7- 44*4 ) 9*fBfgIil* 
*998*«99«»*MMA) -i.—*) 4 1 XX008|g||i* 
**3*8t*0*0M7**4-4)-)2ZAMAM«888i8|i 
♦••§99098 ) )4*44*1 Z7.ZX899891 £?*?••?! 
♦•ill*"* ) 1)4441 4)X4Z7a) ) ) 1 )4llW004 
♦ M 90X)*1 1 ) )S38S«****S47 ) 1 XMA7 ) )M*»4 

******* ********* * ***** ****** ****** 


******************* *************** 

*9*AXA7HA67.A)tt 47X1 AMA777/.71 7X808* 
♦9898MZZ0HX1 1 1 ) + 1 Z) At X W 1 7Z99888* 
*5990*0XXAW17,Zl=7t 1 XM71 ) 170908*** 
*80«0908*Xl XAXX4771 AMA7XXM9988998* 
*0008§0f 807AXX m 4 1 Z)ZX7,X.M000000090* 

♦9809999989^ * ‘ — — — — — 

♦889001 

♦098*9t9999«i9SZMAM9i9i99»99i*9*9* 
♦•00009*8000i00X)Z70f0fif08088888 ♦ 
*99i9999f 999990000090988990999899 ♦ 

*!KSUM**!i*» 00 AXZ» 00 * 099 i!| 9 J|J 
*0999990999000^1 «Z1 7)4X80890919 00* 
*0i!|jJ0fef0MA7.4l a 4 ) ) ) X 000000110 


*«901 


*fiiiii 009 WA)a)*))rf? 7 j« 5 "???lll! 


| 09898 «X)< 


(<0 


*§000099911 D ))4tZlll0»888.) 7.9*080 
♦8000H9 l 71 1 ) ) ) 1 ♦ 1 X ) 27* ) ) 1 1 1 4 ) Z«*9* 
*» 0 f ) S4Z71 1 44 ) ) n ) )«) )Xl7.MWtfY )1 M** 

♦ ************^»iyt ***************** 



Pi 


?.<, « Rocoort of the riot tern , ***' 

(b) OPRA; ®Ha 9 ,in 3 Qn 424,^«3 (c) V \k t { 6 ) OP vA. 





79 


*.^****«?«*t***tt***#**:t* 4 ****** 4 * 


’fffffiftffti 



**«*£) i + ) i z+ } n vn u )-*) 

Jtf *f**A 1 X«»ZXZ=XX 1 L 1 ■*♦ AMPff if »$* 
**i #»i» t# i = ) i z a ) t zx ♦ -« i zmbih ~ ¥ 

***•«! •ISJgzj I *j i 

IffHAZZl )-»lZ*§ 

Jiiz- > X 9 «fit|i»»a«i 

s«js;s«hhh»:iT»;i 

*®iitfi»if *««*** i i — 8Z<*»taai|«t 

!H f 2 I* # 2 *!!K**} ** J >mxAx 5 S#|l®*t* 

1 1 )=ZA -*-+1 *£*##§* 
*i*«e»«BMXZ)*XI — a -= +=+=++! «•»■!* 
*®if#H7,ZAX7.l3---sas4t l)‘/1 ) IZXHlff* 
A) i+=a*-S3t 41ZZXA l«A7t >A«#f * 

s 

***************** ***************+* 


(a) 


****************** **************** 

*«*xx>is)j7.4-mzzi+m«-anA»z§»»t* 
*®t»*»A7Zl lUHlftXI t ! 

*f ® tf »f» A 1 2 -« 7 X 7 = X A 1 74 «4 4 A M !*f»«*®f ® * 
♦••ffflMlAi*) t ZA > 1 7.74 *41 XM»*«® #f f»* 
*OfV9f»»M7. 1 7 ) X 11 Z=* ! ZMMPfff »f ft ft * 
*••••§••<•*•9 A Z 7 1 4 — i 2«®«®f®ft®®fti»* 
*ft«®«®9f*f ®f *8* ) s-1 «•* fffiMf 4P«9* 
*••••• ••••5*»»7-4 ) Af® % lift *••••••• * 

*••••*»•»••■§*»+*****§*§»»»»*#**§§* 

*®V®®®«®|«®V*®M®®®®®®®®® MSfHfi* 

♦fffffffwffff'f fM7X*4 Iff fffff •#§§#♦ 
*if ) i *— *Z»»® if f **»• * * 

*«if»»9B«*e£AX4-1 4-)*l ¥*«f ff 

* 8 ®f a«® 8»8 9 X^)snUl)» A*««»®Sftf® * 
*®»»8®»ft**7 >1=1)1 )=7 A — ~H44«f»®9* 
*»ffffffaZ 7 >aZt-sa -54=4=++l«»iii» 
*ff if »71 A 77 ) a— *sa+ 1 1 1 zl 4 1 1 7 Miff* 
*««B7=1 1 l3-=4=-=l4AZ7XMNAXIAA«»f* 

«• M4. 

********************************** 

(to 


" •• *1 


*<p®*Vf®«S«*a®ft« 9 ft#»* 9 «v *¥*«*¥«** 

*?»#»♦ vvffiifi*** 

♦ W*» '«♦».***##**** = •!“»<- 

♦ 1 •* "♦ a ♦ 7 • • . v ' ~ u V i< . •' • ^ » ins ifi i| t< i v . w t; » 

♦ *fr»l *< " *?<.•' ! ,v 

N >■ ' 44 */ /'• '.epa+B-'.i i '• ,-. ' 

t — ' <tt9al >' *»•• Mj4i v j 1 ,;?• w-a f* “ 

< 1 . / «. 1 A '• H.'d*- ^ J ! 1 .•.* . '■ . ■ «t r 

*=!.., <’*!■, O 1 | 7 ‘\ X i* -,==♦)’ '.1 A ' » 

+ = v . c*4 I * -ft**' 1 , I \ i )4Sr = 1 -.. V \f- ft 4 + 4 j '/ * 

4 r. 1 1 4 4 i Ka» I !*. .■“'I ♦ )1 i V , v * • , ^ 4 , l . ) , , 

4 - * 1 4 * = 7 ) 4 o * 1 1 4 I V < * f. A * ■ 4 + + 21 l' ,1 * 

>*# m *w 

■fc £ * i r: ) s: £ ) \ ♦«^5f i 1 ♦* ) ' j > > * * 

i « * 4 t 4 * f * * 1 * if. M * * * * * * f f * * * * M' * 

(a) 


♦ '/ * ->f »• -7 ® i # 

* *««V*V*f V4«¥»f «*»««¥ Vf49***9^ 

4+av 4 iP ^^v-,»aa.-in; ■>*** ■* ■.«->ft^4<s,fcr .,’• ?nu;* * 

>3 ''>*<» -’if* JH A ft 44*4^ ■ •, ’4i •’ ■,e»S»+‘,. l 1 

* -.X '• », ft W»is:a= "' '^a. : tCi4 4:1:= , V . r* -» = »!•* > 

♦ X"* • *■ .» '*•»' = ■• 'H-#MySd+^' \ *.. 4 -.!sae ■ • .- 

t'S, ' m'Am? . f 'm.-M-Jt v .. v : . l 

9 •’• 1 «. V .• •*,> * •> . V «> . • .' • • .'••■► 

* i ' •»l < f ) '< ,*. ‘ ■ '•* .i ‘ dr = , ■> .' 1 . ,4 r •- V ,i *• -r 

* >'* J 1 < ■ \ A ‘ ’’ ,t A - -i" H • * r » \ .'. ” . ••' t. 

t/,aa1 j v j v a •••.'• .:«i r * + * )'v * 

4 » a 4 1- 1 ) > 771 * 1 ' M » i fl l r ■»■•»-•; > ; t 4 ” ; 1 ! 7 » 

4 } } =+-- f.-ra«--= i 4 = t ■» tt * » 

r ^ 5 ? **** 4 *J 7 *»*»*» ijT.*** \ 

*) S — s-s- ------ It 

♦ f * ♦ 4 » » 4 * * *■ T »• % <•■> **.»*». f» i •( m 4! i- » * » ■( 

0 >) 



Pi/r, r > « 10(a) Hocoa^t ruction of the pattern: i'ff /r l=10,iC2 ! U473, ^"3 • 

(n) PR A, (b) OT'RA • 

Pi". S. 10 (b^ Reconstruction of the ( ?rndefl square : ^JPK «1 1 , ;ILQ* 5 17 , »z 
(n) PRA r (b) OTftA. 
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M * » * * * + **•»♦** » » } * 4 4 * « 

*V**#»#»¥*9*9ti9HV9?H99¥wV»y#¥99SV«' 
*«9w*ifvv»f9Hiiiyu9aitt««*«fv*#y*vu¥^ 
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" M+J 




r*N+J 

>3." 

,1 

r:f N*J- 

• 3 


03900 

04000 

04100 

04200 

04300 

04400 

04500 

04600 

04700 

04800 

04900 

05000 

05100 


6 0 


C 

c 

c 

c 


s,*U Juj-i » n 
T,j= ( T"*l 3 
B(TJ)*.45 
co«TiNUt; 
nu 601=15,18 
00 60 J* 15,1# 

T J = t T - 1 ) * i+o 

TF(B (U).iSO.. 93)8(1 J)=0.0 
I P C 0 ( I 13 . i n .O ,45 ) B C U ) = . 93 
CONTTrlWf 
CALL, 1>RF0CB,M) 

STOP 
F! N D 


r 

C 

c 

c 

c 

c 

c 

c 

c 

c 
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cOTgftTtg# |^R|caHST^yc'DiOH or a .■c&AgjfiO n«ww siuy*. ^ 4 t 

*4 ***f *|f #*#*$♦ *?* **** i* **♦##** ♦***♦&***♦* / **♦*#•* + 

*# ♦♦****♦ *#* ****M********#******f** 4 **♦*;** ♦****♦*» ******¥* 4 * 0 #* 

0 IS fit MATRIX • Wft&fHS TH® PICTURE IS <TK«»5RAt8n 

1 i, D2/>3,04 bM TH® MmXfXBS imn, 

f»JMRMSlO« <m<*24) ’ 1 \ '■ ■ ••'. ' 

■ #««?$*$ ■, . , ■ 

n i *0 , 1 02* • Sf* 103* . t>7 1 04*1 , 00 

lie 

stop ■ T :: 0'-yy-',; : ' : . 

(?(jn ■ ; -," ..;-:;V.' :-,:■ ><■ . ... 

’ “ :y® SQUARE dfittfiRATES' fHf SOUAS®, 

SUHR0UT1 >50^1(6/^,01,02,03,04) 

0 *$l|Pl M 024 : ' 

\miuii ». ,, 

1 0 ' 20 ' ■ 

16)60 TO 10 

§|f4o ; tb/:.4o ; 


m 
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^H^^')!i! *V, i’rCif*'*! aJ «/< k» <«W X >JH(,i‘„ uF" fMfe. AhRA? t< 

J'JW / jW f «‘?i i< 1 1 *- Kf * “iaXItOI WM l,K\/tubi> h'Uu 

or u VMF < 11111 ), 


*) i CU^JIr. k 
ALLtDjFjj t p ti* 


, , , * t Ai’rtX < 1'<J rtft, 

islIrtKUtf ri’Ki* FICOU J‘ In, rt, ihm) 

DlMfcuSIUil bU<>24) 

SUMtlln&l Jw L'H 1 3 U, /) , 0M(2i } 

OirftNbJL'J 4 HtiKCHU ) 

IU 1‘fc.GFR GUlf n,Uu«U \ ) 

‘ j ata Ga / u / 3ua,) f ot3 ( t).G(4a) / o(b,n,ut6 r ij,Gt/ f n 

2^(17 f a , u ■ ( i h. n , oca f i > , « C20. i ) ,(.c 2i f t 5 /h , 1 1- , i h*. 
1 1 ^|y/ ^ ^ ^ ^ ^ r ^ b <1 J f I HO 9 1 <i. \ f 1 Hu § i HU # K H 1 # 1 

^ ^ J ? 4 ^ C * ^ S i i ^ ^ C ^ * f ^ ^ \ ^ I ^ t J * t l * i $ & t 'd ) f is i l ) ? S ) # ^ ( i > f 

2G(17#2) * A) f^Ciy#2i #0(20#^) #ti(2t * 2 J/i i , Ji -i f i a , 

3 1 H * i H tin ftd r i i #■ 1*1*** lH^ f ih*# l*1f f 1 H+ r IHt , jtaX f 1 a 4 9 


41 HX/ 


DAtA «($ t 1},GU # J) f C*U f $>,*(4,3) „ 

^a? ? 35 # GUft f i),v;U9 f 5),a(iu. art' an ; 
3JM » iH , IK , in ,lrt art fli< fill ,lH,,ii*=,li,)iti.,lH., 


"(5, i) ,f (*->,3 1 , UM) 
*» a 3 / 0,1 v j * f i } f i > i 


41H./ 


i&lH, i if 

n,(,cio,D 
iut-a u / 

i s <1 ) 1 " if 

flsktfi) f 

/) f .a >,a 

i ' a , 

1 A f i f 

r'-Ufi), 

■3 ) r i»t l > , i ) 

a<» a-* 

1 I. f 1 rt 


» 

« f 


» i’ C 8 f ♦ J f 
4 J f a I o , t ) , 

i * ' f i i f 
l fi # l •> it f 

f tj l h f ■> J ( 

bJ ,v,t Lo,S> , 
1 *•* f A n , 
liibf l»(^f l H* 

) (» 1 0 ) <) ) f 

t>3 folibfbj j 
1 ii fit' f 

it) , 

f 0 1 4 r / J f 
/Jf&cib,n ( 
t’i flH , 

!•* fit , 


M I (I I f 

OATA aU,4>,GUaJ f b(| r 4J f G(4,4),G(b,4) f G(b # 4J { vj( 1 , 4) 

2i#{i7|4) # Ct C r C 1 0 # 4 } f Gl?\i,4) ,(, 121 .4)/iH a«i a*» , 

31H ,IH ,1H , H an flrt , 1 H , iu ,14 a« ,iit=,lM-,i -(H, 

41HH/ 

1>ATA GCi.b) fGU,’>J fUtifb) ,(,(4.5) ftcC b,b) ,G(.A,b ) ,u( /, i) 
lG(9f5),G( ju.blfGUl.b) ,uU2.b) .6(13.5) .0Ua5) ,(,(13, 
2Gh1|5S,GU8fb5,^tH b ), (,(20,5), GC2{. 55/14 ,1H ,1M 5 

3 in , irt , 1 H , i rt , 1 ft ,irf , l H ,1H ,lti ,|h ,in ,1(1 , 1 nC , 

DATA G(l,t»3#G(2,t>),G(i f h),G(4,6),G(5,'>),r,«j,D).v,(7,b) 

1GC9, 6),G£ 10, o),G( a t, b 5, (,112, ■»), GU (14, o),G(ib, 

2GU 7,t>),aa 8,0 J ,-,(14,0) ,v,(2u.t>) ,U(/l.b)/lH ,ih ,1H , 

3ih ,1h an a -t a« a.t an a » at , ih ,ih a«i an , 

41rtV , irtV/ 

DATA G(i,7),G(2,/),C^(i. / ),(j( 4,7), ,(b, / ),(, (t),?J,v,(/,/) 
lv,(4, / j.GUo, /),(,( 1 1, ?),(,( 12,7). ail, 7 J ,<H1 1,/),G( 15, 
2G(n,7),GUrt, /),v,U a a ,G(2J, 7} ,(,('/!, ;j/lt( a,i a'» , 

31 H , 1 tf , 1 <( fit! flH , 1 rt flrt flrt , t H ,aH ,lt fit ,1*1 , 

4 IMA/ 

DATA D£h/.y, ab,, 22 ,,?b..29.,3i i.3 / 1 15,.b 3, ,5o,U,b, 

1 u.4,.>»/,.7y,y > ibfjbVf,9Jf,47fl,v/ 

DATA NDM/lflfl rl f 1 a f if 1 a f i,2f2f2,3,3,4,4,b,bfOf // 

DATA CHAR/0.6/ 

DATA IttDArtK, UCAH/1H ,1H*/ 

L11KS0 

IP (LlrtF. IF.. ))WI a bo 
IF(Dl^h..£s)« <> ) t»»( i j i t i 
MDI«fe.= UD 

XF(LIrtPS.£v u ) "tDi <S=F0 
XrtIN=B(U ?X«AK=srt(U 
xsuw=o. 

DO JbJal, v 
DO 160*1,4 
1J*U-D*«+J 
XSUM*XSUM+D(IJ) 

IP CXMlf»-S(U))Uf A2,i*l 

iFUNAx-btu))ua<>ab 

XrtAX«rt( Ivl) 

CO 4TX A*Hfe. 

ht*FU)AlCrt)*CMAK 

XDlSfsP DO U'Ui) /FLOAT (Wi + l) 

iF((x«AX-x«iiOfG(a.p:-io)GO ru 18 

lFCDI*fi..rtt.D)F«l,a 72 

RcTURM 

IP ( C^ + 2) .Gl'.bLl lc.)RP tUrtrt 
FKAC=l,/(X3AX-XrtX O 



If O.G'i.lUuJGu T«i 72 

NON=tbO- . I J/4-1 

IKCLdiP,. Rv», 
nu ?oi=I,.i<i* 

PRINT bh 
CUimwuh 
DU 241*1 , J Ju 
uu 24J*1 { 1 
UNCI, JlaJuuftWK 
CONTINUE 
X<4IN={iD,l „ fc - ,i > / ^ 

I*AX*l.M4 + w-l 

II1*I«I*-1 

II2=lNAx+i 

DO 261*111,112 

LN(I#l)*Aa»AH 

CONTINUE 

PRINT 74, (MCI ,1 1,1 = 1 , 

PK I w £ flb 
Y»PLUAJ C'i) + .5 
DO 46K*1.«1 
ysY-XDlbT 

ji*y 

yi=flua'i cju 

J2»J1+1 

11*0 

LMAX*0 

DU 401*1 IX'i,J.4AX 

11*11+1 

iui«ai-u*^jt 

iu2*ai-n»^j2 

D*c0auiUibUiU2)-<sciium*ty-HH num-v t 
DO 2SM*t ,2i 
IFCD-DENC 4>)*u,4b,2P 
CONTINUE 

T*(DENCH2 ) +ORN (nl 13/2. 

IFC0-T332, 32,34 
li»Ml 

GO TO 3H 

L*M2 

GO ft) 3H 
L=M 

DU 400 = 1 , 7 

IF ( ti.GT . M* AX ) h iAX=L 

LNC I , J)=GU< , J) 

CONTINUE 
NJ*NLN(L^AXJ-1 
if CNJ.Lfc ,0) Gtl 10 44 
DU 42 J*1 , N J 

PRINT 7b, (D id ,J) ,1*1, Mot,) 

CU.*1 !• lM 

NJ1= <U + i 

PRINT /cs, (DnCX , N J 1 ) ,1*1, i'j ] K) 

CUM INUfc 
DO 481*111,112 
LNCI, l)*T r ' -R 
CONTINUE 

PfPT 74, CIA U ,11 ,1*1 , ■ * I «e ) 

MI)f,Ijr*NLI 4b /ll 
flUEG*22/i!l(Ku 
DU 561 1= 1 1 , NtiFG 
DU 5QJ = 1 , 7 
00 501 = 1,1 IN 

UNC1,J)=IAuANk 

CONTINUE 

DO 52J*i ,11 

If>UBl = (J*‘l)*lt+t 

16U&2*Cu-l) + Ul~i)*'<LRG 

IMJB3=ISUH2-*i 

IFC I5UD2 ,DU»0 J J bUB2*t 

RLfcGCJ)*?f)fcMUuD21+DE 1 *abUb3)}/(FPAC*2.)+XMlN 
DO 52K*t , t 


)*FMC 


i^jQfauai ,-o=Ju,sa.u,/u 

cumimuk 

W) 5U*l,t> 

PhI«.T 7o r (hw(J, U, 1=1, 
PRJLmC 78, u C l, 

P*I*T rtO,(RLfe.aCJ),J=l , 
continue 

it (lire. ue.o)«j si 6i 
DO 50J«1 , 7 
00 b8I=l, 130 
LU ( X ,«J) = IRliAwE 

cum Ihuk 

RLE(,(1 )=0KN( il l/FRACfX' 

DU 60K=l ,7 

ON C 1 , K ) =G (21 , K.) 

continue 

PRINT 68 

DU 621*1,0 

PRINT 76, It 1(1,1) 

PRINT 78,64(1,7) 

PRINT 80,RLe.fHl) 
CONTINUE 

IFCLlUK.UE.OlPfti If U 
RETURN 

FORMAT ( OKI* JVF*) 

FORMAT (/) 

FORflATCSrft* W4*} 

FORMAT UH ,30A,10£)a1) 
FURMA|C1H+,3JK,Iu«»a 1) 
FORMA! t 1 H t 3QK f IOOmI ) 
FURMAT(IUU,Elo» U) 
FORMATl/) 

END 


• ILL ■>$. 
‘ihl .h 
•LEU) 


«ll * 


I © 

20 

40 


60 

70 


80 


100 

fto 

UQ 


5 UhIoIJTT ~ f :’ 7 5 ) h“jh u ' 


KLFKHfc .- JCtr.oS 


a i f*t» v t l ho i . i jcl r«.r p'liju * Ji 1 " •>.).'* 

C f‘V. PTC J Ot-t ' A 1 H 1 X . 

F f‘»fi *-» rC :» ,MM CTS.U I'iC » OHfc 

U i* 4 l -> fV - L 1 KI P * I Hr . 

b t’ Ifc WfC 1 »J fJ MCti <M PS Inrttf 

d i‘ C T 1 r ' * 

y “i tu 1 A A P,Lr ftjt THb F.vuJfCJJd 


# ' * , 

til f I p '0- 
l.S ti \ l . U ’ 


S0BRUU1 CAbbCdi 

PROJT t SHHPrtU.il I f 


I 


red o * 

.)£$} **lo A Gi-b. 

P«if»y ; Udf’AlJS l * Kb PRuo C i *’i 1 

QOTM ; PKHF HR''.'. Ui* t -I'l'Jlt ’^LlZA H'i 

FtC'Hir : PRl*7h u »l l 1 'I II A«,j* . 
PRH ! CALC* LAI Kit "Hud* < HM »ft, 


* 1 


f . 


> i 


/-■» * , U >.' ') L * <• !') 


SU6H0UTIAK PRrUCC, ) 

DIMENSIO I Aao24),MAU0/4) ,*■ C * 

PKAL C(1024),n r Art f iFflu^A) 

INTEGER P , PH f & 

REAL LC20) 

PRINT 1 

format? ir , 'KfeCU'-i&ri'Ucri no uw'itftiWt^u.btun. .7/) 

PI*3,14?8570 

MPR=§ 

DU 5LL=1 , 4 
MARK*LL 
TYPE*, MARK 
IF(Ll««LT.2)&0 1U Id 

LUjsiU JLm = Q .? LCn = 0 .; L ( 4 l = o.;l (*>)=( . i / i = 

Lc5)*o,;btlO)*d.;ul U )=(,,?;>=<» 

THETA*©, ?5*l 
IFCLCSJ.l.T.O,^} ,1t. X i 
IKLCS),bl,3?,J'. J ?n 
CJU T«1 7 U 
LCi>)=b(S) + .7 

CALL PK‘1Ji ( I . I K Y 4 , C ,!<»,<, A, B,K.,P,5,1 
IF(P,LL,1)GU 'Id 20 
IF (MARK.Ey, 1 ) ”‘=‘ *-i 

DU 401=1,10X4 
AA ( I ) =A ( 1 ) 

HB = H 
PP=P-1 
DO 601=1, PR 
KA ( 1 ) = K ( I } 

THfcirAsrd'lH+PT/O. J3=S + 1 
IF (S.GT, 9 )GU rn I 60 
IFCL?S) ,GK.32)GJ TO /o 
L(S)=L(S)+,6 
DO 90I=t ,10X4 
A C I ) =0 , 


.1 n = 


EalPrKOJHI i * ,b,K,P,S) 

iEU ' CSJ ^ TiJ / j . >< JO 7 o 

1-1 


IFCP , LF ^ f)Su Til 00 


)/Z*AA(KF(D) 


CAijL Pi'LOUA,F f y J 
CALL LX hA.AA.Z) 

iFCz.feo.o.iLO to no 

DO 100 1 = 1 , PR 

f Ci < KU ))= MAKcn )- n-f . _ 
XFCF ( K^t n ) .0 I .< ,)Ft K * l I J1 * 0 , 

fit J 1201 = 1 , 1 02 , 

AA( I )=A( r ) 

HH=L 
PR=R»1 



140 

lbO 

1/0 


80 

00 


220 


240 

300 


6 

7 

222 


244 

5 


DU 1401=1 t Pk> 

RK(1)=K( I ) 

IiUS-0 

i/0 i / 0 & J f ^ 

iWr* S,uu 1J ?u 

CHUTI/Uit 

FORUAUlH ,i 2*5,2) 

TUFTA=0.)S=1 

LC8)=0. 

LCS)=L(S)+.5 

CAUL PR0JTC IHhlA »C,ii, «,A,R,*,p,S) 
IFCP.LE.UGO TO 220 ' ' ' 

CALL QRTH(aa,BB,A,B) 

CALL PROD(AA,F,$) 

CALL PKUDCAA,AA,Z) 

DO 240T=1,PP 

PCKKC pJ=FCKK(i)i-C\-Lrt}/4*AA(A} ( L)) 
XFO* CKK( U).Lr.O. )t (Kr CTJ )=u 
PR IN f 300, MARK, NPR,Hfc.U 
tOKMAI'CH t *F‘)R *,l4, ¥ i t'kHA 


p ufinn* uo i r 'in 

I'tyUATIUuS V) 

CALL PICDUJ'CF, 0,1024) 

MISsO 

DO 61=1,1024 

1fIf(?).L1..30)^FCU=0. 

if Ml :»::i85J#?iK!is b,t 

00 71*1,1024 

IF(CCI).NE.MFCI))aS*aA+l 
PRINT 222, MIS 
TYPE 222, MIS 
FORMAT UH , NO, 

CALL ERRCC,F,<0 
PRINT 244 
FORMAT! 1H ,///) 

CONTINUE 
RETURN 
END 


» 14 , ) TLRALlti *lb) # , 14 , 'PROJECTIONS 


,&«)> I ( 1 )=, 48 


)F • UbMAlCHtb * ' ,15) 


9 I 4 f 


utm 1 
W r 

t»07.)' 

Hk>8</ 
“lit- 

l4t,r 
.5 J(' 
„ 6 ( ,f ’ 

1171, f 

020u?t 

Q2lon 

O2?0 r " 


P p f^P^#Pw7ou**** 

njwpoyoM"p?a024nRR<rjJ4r 

1 ‘« 8o«Joi?pf l??RMI> 

®iTUR#» 

r— 

c * mmnm ~i 35, S5f hosSJ AEm^Qt^""’" ,, ’ 

' “y«m«*iwi 

M®®’ 

jo AAfn»AA(n*irz/zy»iviT) 

R8*«!f|-vZ/ZY*a 

giJ^PN 



M 


Ou 1 v 

ou*> 

OUT 

oik- 

QO*" 

OOf' 

00 7 
Oo lJ 

00 ° 
OK ' 1 
011 < ' 
:oi'a • 
onuf 
01AU 
01«5* * 

oim 

017 » ( 
OlPv '* 
OiA o 
02 <K'> 
02*00 
022 UP 
02100 
02400 
O?50O 
02600 
027 W 0 
02600 
02*00 
OK 
03 ] 

si 

o! 
o; 
oi 


r>, 

r~. 



04A^O 
Ml ' 

«» i< 

0-*7 - 
04 f<» 

04 A < 
Ot »0 ' 

Obi v,< 
Ob? ' 
05M 
Oi)4» ' 
Ob*>h< 
Obfi r 
0&7 o 
05 fU/f 
ObA ' 
ObP 
Obi ' 1 
06?>< 
Ob3< 0 
064 v p 
lto5o? 
6600 
>|700 
>6000 
»00(< 
070 of 
0710^ 
§7200 
07 lOO 
07400 
07500 


ffOBrCnUTtHfi i»RACC • »> 

ni^ 4 ),FU 0 ? 4 ),K( 6 OO),fitlO 2 *),MVUU 24 l 

T m* & h f r*lr#5 
Riftu Ltse$,fL 
w tm*6 - 
* *3 
PI* 

no 

f A§ 

*<* 


1 0 

20 


1 .< 

1 1 

1 / 


f 

7 

72 ° 

5 
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*O 7 tC 3 )* 0 ?fcf 4 l*O?r,( 5 )* 0 !f£ J rb)*OrL( 7 )s 0 Ju(e)*O 

J 4 sulS *0 

* l 5)€U 'to ?y 

* ‘ ' «0 


A>P»K*Pf 5 l 


i)WEy* *f .)♦! 


0 70 
as 1 *Pp 

-f Y-i VZ*Mr(T)l 
#>F CKfI 3 )*U 
,/b. ? Sa* + 1 
GO in Ib0 
,32)no I*n 7 *• 

tf v 24 


pftQ.Jf(THtTrt,u ,u, ,*, 0 , 1 ^°, S' 
# 01 * 32)60 7 t 
« K :.<ko ir pu 
. WRK^v.Hi Fg*’^\ 2 i 
pp*p*i < 

pmu pRtn( 4 ,F,Y) 

” a I L t'Ryn ( ii 1 A| Z) 

TMzjtn.ojno tn i«o 

i, 11 UTat # Pp 

r v r jlsFM (Tn-Q<**KVZ*AiK(T>) 
T fO(ftfi’),r.p,)KK(i)jso, 

r t '’t'l 

’ i 7 1 "SSI ,f 

K <- ,r.') .ne *i 2 )T,SSsT.*s+l 
n* K^b.t M>ni 10 it«o 


K # 4 )< 

ri n 1 „ 



P| T 

i * 1 u'> 

1 ^ r 

P» T 

p 

7 n 

F ( ,4 

if 1 f i W 


f\ r 

. M ( K’ 

/* M 

r a 

A if//) 

r 

^ *T 

a PT C* 1 <- 


Pul 



r Jf« 

55 ^ 


< j 

t * 


tP 7 0 


F'M m*?K 
i Pk * M 4 , 
1 5 . 2 i 

. >’ , t «,? n 


NFO * MS. 


ITgRfcTinw(S) * M*/„) 


Pi OTsl^t/K 
* r M*'> 

is n f j i,r b ., 5 j“mj*i 

^ J H's* , U>< + 

FlT)li.l 6 «hI 34 | 

Pnf M ?,2,U« 

M- i ( i 

r.jp * " M * , 1 
t"V T jT 4 yO 

Mia kPK(C,t »«0 
ru n l'TN”L 
PF^JR^ 

P* *0 


‘*U. rj> ^3 1 ATCHE 8 * M 57 


f)u 1 r - 

oo n 

Oo* r- 

Ovi • / 

oo r r 

0)f' / 

007 , » 

' 

OOQ / “* 

010,7 

Oil*' ’ 

on, iu 

on,> 

01 1 j 

cn*j ' 

o lfi. ?o 

0i7^ 

01" 
or» « rt 
0201,* 

021 V * ' 

02?o0 
02 3 or 


02400 300 
02*>u0 400 
026yO 400 



0‘i 7 
O'* 1 
0<*S 
0«7 

o<n> ■ >r * 

0*0 f . 

OW Rb 

Obi 
05?> ' 

Ob ?■, 


C,, fc^iSfC, 



hrrrrw 


71 "ft. '-ST DM r(10?4) ,1 
VAi. 1R&RR 
' >,*& * n 

c ,u"0*C.;^J U=0.;<5u Vs^.ylUMjs^ 

C0?lTNU J e 7tAliSCCn) * rfl3f 

AVG»$UMO/iJU 
ou 2oi»i.v' 

3UM4*5'NH(CCn+flVG}**2 

ot,r,T&#sofnrsoM/%o 

'.7rt{Pa3uM2/S'T ,0 
VAOsSuMJ/ , J 

vsasosrrs’i.n/su”^) 

PRTtfT loo 
PmTr? 4u0 

°RImT 5wO,l>e,fjM,*fREHR,VAR,P'ir 

FORMAXQH , * THE ERROR mEaSURE* ARE I *,/) 
FOPrtArUH ,*011X4 ■ ,F8.5 , * MR ERR * * »FR, 

t RHS a *,F8,5) 
f 1,1024 

).G£..5JFF(n»7 .0 

f I7a0, 

SUHHO«0, »30H M lsO,,-5n M «2sO. r$n«¥jsO # j su*«*4«n, 

tj<4I«?uA+cu) 

uAfAUsu*»'.i + (cm- B 'i m)*i? 

■‘*2#SUM«? + AR5(' , (T J-PF(T)) 

1*T HE 
r,«SoMfn)/" ( 
f)0 25l»l,r,' 

4« *H4s4il' l+fan-AKj*»? 

<?L '•* IsSbM i7 + f «•*’( I )-AV(,)**2 
\. r ’b‘ n .\a.b?rt M v4ir^w1 /fl„l 

* *ERR*4RMu2/bH < 10 

* )R r(S>» 1-tl /SU*M4) 
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